限制对比度自适应直方图均衡化算法原理、实现及效果

Matlab 同时被 2 个专栏收录
34 篇文章 0 订阅
59 篇文章 0 订阅

http://blog.csdn.net/laviewpbt/article/details/8769877#comments

一、自适应直方图均衡化(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的斜度因此,也限制了变换函数的斜度。直方图被裁剪的值,也就是所谓的裁剪限幅,取决于直方图的分布因此也取决于领域大小的取值。

     通常,直接忽略掉那些超出直方图裁剪限幅的部分是不好的,而应该将这些裁剪掉的部分均匀的分布到直方图的其他部分。如下图所示。

 

        Clahe-redist.svg
 
    这个重分布的过程可能会导致那些倍裁剪掉的部分由重新超过了裁剪值(如上图的绿色部分所示)。如果这不是所希望的,可以不带使用重复不的过程指导这个超出的部分已经变得微不足道了。

     2. 通过插值加快计算速度

        如上所述的直接的自适应直方图,不管是否带有对比度限制,都需要对图像中的每个像素计算器领域直方图以及对应的变换函数,这使得算法及其耗时。

        而插值使得上述算法效率上有极大的提升,并且质量上没有下降。首先,将图像均匀分成等份矩形大小,如下图的右侧部分所示(8行8列64个块是常用的选择)。然后计算个块的直方图、CDF以及对应的变换函数。这个变换函数对于块的中心像素(下图左侧部分的黑色小方块)是完全符合原始定义的。而其他的像素通过哪些于其临近的四个块的变换函数插值获取。位于图中蓝色阴影部分的像素采用双线性查插值,而位于便于边缘的(绿色阴影)部分采用线性插值,角点处(红色阴影处)直接使用块所在的变换函数。

        Clahe-tileinterpol.svg
 
     这样的过程极大的降低了变换函数需要计算的次数,只是增加了一些双线性插值的计算量。

     CLAHE算法的源代码参考:

[cpp]  view plain copy
  1. View Code   
  2.   
  3. /* 
  4.  * ANSI C code from the article 
  5.  * "Contrast Limited Adaptive Histogram Equalization" 
  6.  * by Karel Zuiderveld, karel@cv.ruu.nl 
  7.  * in "Graphics Gems IV", Academic Press, 1994 
  8.  * 
  9.  * 
  10.  *  These functions implement Contrast Limited Adaptive Histogram Equalization. 
  11.  *  The main routine (CLAHE) expects an input image that is stored contiguously in 
  12.  *  memory;  the CLAHE output image overwrites the original input image and has the 
  13.  *  same minimum and maximum values (which must be provided by the user). 
  14.  *  This implementation assumes that the X- and Y image resolutions are an integer 
  15.  *  multiple of the X- and Y sizes of the contextual regions. A check on various other 
  16.  *  error conditions is performed. 
  17.  * 
  18.  *  #define the symbol BYTE_IMAGE to make this implementation suitable for 
  19.  *  8-bit images. The maximum number of contextual regions can be redefined 
  20.  *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256 
  21.  *  contextual regions is not recommended. 
  22.  * 
  23.  *  The code is ANSI-C and is also C++ compliant. 
  24.  * 
  25.  *  Author: Karel Zuiderveld, Computer Vision Research Group, 
  26.  *         Utrecht, The Netherlands (karel@cv.ruu.nl) 
  27.  */  
  28.   
  29. #ifdef BYTE_IMAGE  
  30. typedef unsigned char kz_pixel_t;     /* for 8 bit-per-pixel images */  
  31. #define uiNR_OF_GREY (256)  
  32. #else  
  33. typedef unsigned short kz_pixel_t;     /* for 12 bit-per-pixel images (default) */  
  34. # define uiNR_OF_GREY (4096)  
  35. #endif  
  36.   
  37. /******** Prototype of CLAHE function. Put this in a separate include file. *****/  
  38. int CLAHE(kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes, kz_pixel_t Min,  
  39.       kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,  
  40.       unsigned int uiNrBins, float fCliplimit);  
  41.   
  42. /*********************** Local prototypes ************************/  
  43. static void ClipHistogram (unsigned long*, unsigned int, unsigned long);  
  44. static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,  
  45.         unsigned long*, unsigned int, kz_pixel_t*);  
  46. static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,  
  47.            unsigned int, unsigned long);  
  48. static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);  
  49. static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,  
  50.     unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);  
  51.   
  52. /**************     Start of actual code **************/  
  53. #include <stdlib.h>             /* To get prototypes of malloc() and free() */  
  54.   
  55. const unsigned int uiMAX_REG_X = 16;      /* max. # contextual regions in x-direction */  
  56. const unsigned int uiMAX_REG_Y = 16;      /* max. # contextual regions in y-direction */  
  57.   
  58.   
  59.   
  60. /************************** main function CLAHE ******************/  
  61. int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,  
  62.      kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,  
  63.           unsigned int uiNrBins, float fCliplimit)  
  64. /*   pImage - Pointer to the input/output image 
  65.  *   uiXRes - Image resolution in the X direction 
  66.  *   uiYRes - Image resolution in the Y direction 
  67.  *   Min - Minimum greyvalue of input image (also becomes minimum of output image) 
  68.  *   Max - Maximum greyvalue of input image (also becomes maximum of output image) 
  69.  *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X) 
  70.  *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y) 
  71.  *   uiNrBins - Number of greybins for histogram ("dynamic range") 
  72.  *   float fCliplimit - Normalized cliplimit (higher values give more contrast) 
  73.  * The number of "effective" greylevels in the output image is set by uiNrBins; selecting 
  74.  * a small value (eg. 128) speeds up processing and still produce an output image of 
  75.  * good quality. The output image will have the same minimum and maximum value as the input 
  76.  * image. A clip limit smaller than 1 results in standard (non-contrast limited) AHE. 
  77.  */  
  78. {  
  79.     unsigned int uiX, uiY;          /* counters */  
  80.     unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */  
  81.     unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */  
  82.     unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */  
  83.     kz_pixel_t* pImPointer;           /* pointer to image */  
  84.     kz_pixel_t aLUT[uiNR_OF_GREY];        /* lookup table used for scaling of input image */  
  85.     unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/  
  86.     unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */  
  87.   
  88.     if (uiNrX > uiMAX_REG_X) return -1;       /* # of regions x-direction too large */  
  89.     if (uiNrY > uiMAX_REG_Y) return -2;       /* # of regions y-direction too large */  
  90.     if (uiXRes % uiNrX) return -3;      /* x-resolution no multiple of uiNrX */  
  91.     if (uiYRes & uiNrY) return -4;      /* y-resolution no multiple of uiNrY */  
  92.     if (Max >= uiNR_OF_GREY) return -5;       /* maximum too large */  
  93.     if (Min >= Max) return -6;          /* minimum equal or larger than maximum */  
  94.     if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */  
  95.     if (fCliplimit == 1.0) return 0;      /* is OK, immediately returns original image. */  
  96.     if (uiNrBins == 0) uiNrBins = 128;      /* default value when not specified */  
  97.   
  98.     pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);  
  99.     if (pulMapArray == 0) return -8;      /* Not enough memory! (try reducing uiNrBins) */  
  100.   
  101.     uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */  
  102.     ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;  
  103.   
  104.     if(fCliplimit > 0.0) {          /* Calculate actual cliplimit     */  
  105.        ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);  
  106.        ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;  
  107.     }  
  108.     else ulClipLimit = 1UL<<14;          /* Large value, do not clip (AHE) */  
  109.     MakeLut(aLUT, Min, Max, uiNrBins);      /* Make lookup table for mapping of greyvalues */  
  110.     /* Calculate greylevel mappings for each contextual region */  
  111.     for (uiY = 0, pImPointer = pImage; uiY < uiNrY; uiY++) {  
  112.     for (uiX = 0; uiX < uiNrX; uiX++, pImPointer += uiXSize) {  
  113.         pulHist = &pulMapArray[uiNrBins * (uiY * uiNrX + uiX)];  
  114.         MakeHistogram(pImPointer,uiXRes,uiXSize,uiYSize,pulHist,uiNrBins,aLUT);  
  115.         ClipHistogram(pulHist, uiNrBins, ulClipLimit);  
  116.         MapHistogram(pulHist, Min, Max, uiNrBins, ulNrPixels);  
  117.     }  
  118.     pImPointer += (uiYSize - 1) * uiXRes;          /* skip lines, set pointer */  
  119.     }  
  120.   
  121.     /* Interpolate greylevel mappings to get CLAHE image */  
  122.     for (pImPointer = pImage, uiY = 0; uiY <= uiNrY; uiY++) {  
  123.     if (uiY == 0) {                      /* special case: top row */  
  124.         uiSubY = uiYSize >> 1;  uiYU = 0; uiYB = 0;  
  125.     }  
  126.     else {  
  127.         if (uiY == uiNrY) {                  /* special case: bottom row */  
  128.         uiSubY = uiYSize >> 1;    uiYU = uiNrY-1;     uiYB = uiYU;  
  129.         }  
  130.         else {                      /* default values */  
  131.         uiSubY = uiYSize; uiYU = uiY - 1; uiYB = uiYU + 1;  
  132.         }  
  133.     }  
  134.     for (uiX = 0; uiX <= uiNrX; uiX++) {  
  135.         if (uiX == 0) {                  /* special case: left column */  
  136.         uiSubX = uiXSize >> 1; uiXL = 0; uiXR = 0;  
  137.         }  
  138.         else {  
  139.         if (uiX == uiNrX) {              /* special case: right column */  
  140.             uiSubX = uiXSize >> 1;  uiXL = uiNrX - 1; uiXR = uiXL;  
  141.         }  
  142.         else {                      /* default values */  
  143.             uiSubX = uiXSize; uiXL = uiX - 1; uiXR = uiXL + 1;  
  144.         }  
  145.         }  
  146.   
  147.         pulLU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXL)];  
  148.         pulRU = &pulMapArray[uiNrBins * (uiYU * uiNrX + uiXR)];  
  149.         pulLB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXL)];  
  150.         pulRB = &pulMapArray[uiNrBins * (uiYB * uiNrX + uiXR)];  
  151.         Interpolate(pImPointer,uiXRes,pulLU,pulRU,pulLB,pulRB,uiSubX,uiSubY,aLUT);  
  152.         pImPointer += uiSubX;              /* set pointer on next matrix */  
  153.     }  
  154.     pImPointer += (uiSubY - 1) * uiXRes;  
  155.     }  
  156.     free(pulMapArray);                      /* free space for histograms */  
  157.     return 0;                          /* return status OK */  
  158. }  
  159. void ClipHistogram (unsigned long* pulHistogram, unsigned int  
  160.              uiNrGreylevels, unsigned long ulClipLimit)  
  161. /* This function performs clipping of the histogram and redistribution of bins. 
  162.  * The histogram is clipped and the number of excess pixels is counted. Afterwards 
  163.  * the excess pixels are equally redistributed across the whole histogram (providing 
  164.  * the bin count is smaller than the cliplimit). 
  165.  */  
  166. {  
  167.     unsigned long* pulBinPointer, *pulEndPointer, *pulHisto;  
  168.     unsigned long ulNrExcess, ulUpper, ulBinIncr, ulStepSize, i;  
  169.     long lBinExcess;  
  170.   
  171.     ulNrExcess = 0;  pulBinPointer = pulHistogram;  
  172.     for (i = 0; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */  
  173.     lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;  
  174.     if (lBinExcess > 0) ulNrExcess += lBinExcess;      /* excess in current bin */  
  175.     };  
  176.   
  177.     /* Second part: clip histogram and redistribute excess pixels in each bin */  
  178.     ulBinIncr = ulNrExcess / uiNrGreylevels;          /* average binincrement */  
  179.     ulUpper =  ulClipLimit - ulBinIncr;     /* Bins larger than ulUpper set to cliplimit */  
  180.   
  181.     for (i = 0; i < uiNrGreylevels; i++) {  
  182.       if (pulHistogram[i] > ulClipLimit) pulHistogram[i] = ulClipLimit; /* clip bin */  
  183.       else {  
  184.       if (pulHistogram[i] > ulUpper) {        /* high bin count */  
  185.           ulNrExcess -= pulHistogram[i] - ulUpper; pulHistogram[i]=ulClipLimit;  
  186.       }  
  187.       else {                    /* low bin count */  
  188.           ulNrExcess -= ulBinIncr; pulHistogram[i] += ulBinIncr;  
  189.       }  
  190.        }  
  191.     }  
  192.   
  193.     while (ulNrExcess) {   /* Redistribute remaining excess  */  
  194.     pulEndPointer = &pulHistogram[uiNrGreylevels]; pulHisto = pulHistogram;  
  195.   
  196.     while (ulNrExcess && pulHisto < pulEndPointer) {  
  197.         ulStepSize = uiNrGreylevels / ulNrExcess;  
  198.         if (ulStepSize < 1) ulStepSize = 1;          /* stepsize at least 1 */  
  199.         for (pulBinPointer=pulHisto; pulBinPointer < pulEndPointer && ulNrExcess;  
  200.          pulBinPointer += ulStepSize) {  
  201.         if (*pulBinPointer < ulClipLimit) {  
  202.             (*pulBinPointer)++;     ulNrExcess--;      /* reduce excess */  
  203.         }  
  204.         }  
  205.         pulHisto++;          /* restart redistributing on other bin location */  
  206.     }  
  207.     }  
  208. }  
  209. void MakeHistogram (kz_pixel_t* pImage, unsigned int uiXRes,  
  210.         unsigned int uiSizeX, unsigned int uiSizeY,  
  211.         unsigned long* pulHistogram,  
  212.         unsigned int uiNrGreylevels, kz_pixel_t* pLookupTable)  
  213. /* This function classifies the greylevels present in the array image into 
  214.  * a greylevel histogram. The pLookupTable specifies the relationship 
  215.  * between the greyvalue of the pixel (typically between 0 and 4095) and 
  216.  * the corresponding bin in the histogram (usually containing only 128 bins). 
  217.  */  
  218. {  
  219.     kz_pixel_t* pImagePointer;  
  220.     unsigned int i;  
  221.   
  222.     for (i = 0; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */  
  223.   
  224.     for (i = 0; i < uiSizeY; i++) {  
  225.     pImagePointer = &pImage[uiSizeX];  
  226.     while (pImage < pImagePointer) pulHistogram[pLookupTable[*pImage++]]++;  
  227.     pImagePointer += uiXRes;  
  228.     pImage = &pImagePointer[-uiSizeX];  
  229.     }  
  230. }  
  231.   
  232. void MapHistogram (unsigned long* pulHistogram, kz_pixel_t Min, kz_pixel_t Max,  
  233.            unsigned int uiNrGreylevels, unsigned long ulNrOfPixels)  
  234. /* This function calculates the equalized lookup table (mapping) by 
  235.  * cumulating the input histogram. Note: lookup table is rescaled in range [Min..Max]. 
  236.  */  
  237. {  
  238.     unsigned int i;  unsigned long ulSum = 0;  
  239.     const float fScale = ((float)(Max - Min)) / ulNrOfPixels;  
  240.     const unsigned long ulMin = (unsigned long) Min;  
  241.   
  242.     for (i = 0; i < uiNrGreylevels; i++) {  
  243.     ulSum += pulHistogram[i]; pulHistogram[i]=(unsigned long)(ulMin+ulSum*fScale);  
  244.     if (pulHistogram[i] > Max) pulHistogram[i] = Max;  
  245.     }  
  246. }  
  247.   
  248. void MakeLut (kz_pixel_t * pLUT, kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrBins)  
  249. /* To speed up histogram clipping, the input image [Min,Max] is scaled down to 
  250.  * [0,uiNrBins-1]. This function calculates the LUT. 
  251.  */  
  252. {  
  253.     int i;  
  254.     const kz_pixel_t BinSize = (kz_pixel_t) (1 + (Max - Min) / uiNrBins);  
  255.   
  256.     for (i = Min; i <= Max; i++)  pLUT[i] = (i - Min) / BinSize;  
  257. }  
  258.   
  259. void Interpolate (kz_pixel_t * pImage, int uiXRes, unsigned long * pulMapLU,  
  260.      unsigned long * pulMapRU, unsigned long * pulMapLB,  unsigned long * pulMapRB,  
  261.      unsigned int uiXSize, unsigned int uiYSize, kz_pixel_t * pLUT)  
  262. /* pImage      - pointer to input/output image 
  263.  * uiXRes      - resolution of image in x-direction 
  264.  * pulMap*     - mappings of greylevels from histograms 
  265.  * uiXSize     - uiXSize of image submatrix 
  266.  * uiYSize     - uiYSize of image submatrix 
  267.  * pLUT           - lookup table containing mapping greyvalues to bins 
  268.  * This function calculates the new greylevel assignments of pixels within a submatrix 
  269.  * of the image with size uiXSize and uiYSize. This is done by a bilinear interpolation 
  270.  * between four different mappings in order to eliminate boundary artifacts. 
  271.  * It uses a division; since division is often an expensive operation, I added code to 
  272.  * perform a logical shift instead when feasible. 
  273.  */  
  274. {  
  275.     const unsigned int uiIncr = uiXRes-uiXSize; /* Pointer increment after processing row */  
  276.     kz_pixel_t GreyValue; unsigned int uiNum = uiXSize*uiYSize; /* Normalization factor */  
  277.   
  278.     unsigned int uiXCoef, uiYCoef, uiXInvCoef, uiYInvCoef, uiShift = 0;  
  279.   
  280.     if (uiNum & (uiNum - 1))   /* If uiNum is not a power of two, use division */  
  281.     for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;  
  282.      uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {  
  283.     for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;  
  284.          uiXCoef++, uiXInvCoef--) {  
  285.         GreyValue = pLUT[*pImage];           /* get histogram bin value */  
  286.         *pImage++ = (kz_pixel_t ) ((uiYInvCoef * (uiXInvCoef*pulMapLU[GreyValue]  
  287.                       + uiXCoef * pulMapRU[GreyValue])  
  288.                 + uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]  
  289.                       + uiXCoef * pulMapRB[GreyValue])) / uiNum);  
  290.     }  
  291.     }  
  292.     else {               /* avoid the division and use a right shift instead */  
  293.     while (uiNum >>= 1) uiShift++;           /* Calculate 2log of uiNum */  
  294.     for (uiYCoef = 0, uiYInvCoef = uiYSize; uiYCoef < uiYSize;  
  295.          uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {  
  296.          for (uiXCoef = 0, uiXInvCoef = uiXSize; uiXCoef < uiXSize;  
  297.            uiXCoef++, uiXInvCoef--) {  
  298.            GreyValue = pLUT[*pImage];      /* get histogram bin value */  
  299.            *pImage++ = (kz_pixel_t)((uiYInvCoef* (uiXInvCoef * pulMapLU[GreyValue]  
  300.                       + uiXCoef * pulMapRU[GreyValue])  
  301.                 + uiYCoef * (uiXInvCoef * pulMapLB[GreyValue]  
  302.                       + uiXCoef * pulMapRB[GreyValue])) >> uiShift);  
  303.         }  
  304.     }  
  305.     }  
  306. }  


 

     上面的代码中对于各块之间的插值部分的编码技巧很值得学习和参考。

      以上描述均翻译自:http://en.wikipedia.org/wiki/CLAHE#Contrast_Limited_AHE

      Karel Zuiderveld提供的代码:

[cpp]  view plain copy
  1. if (pulHistogram[i] > ulUpper)      
  2.             {       /* high bin count */     
  3.                 ulNrExcess -= (pulHistogram[i] - ulUpper); pulHistogram[i]=ulClipLimit;     
  4.             }     

     应该修正为:

[cpp]  view plain copy
  1. if (pulHistogram[i] > ulUpper)      
  2.             {       /* high bin count */     
  3.                 ulNrExcess -= (ulClipLimit -pulHistogram[i]); pulHistogram[i]=ulClipLimit;     
  4.             }     

      同时,各位也可以参考下matlab的adapthisteq.m文件,该文件的代码基本是严格按照 Karel Zuiderveld作者的原文写的,贴出如下:

 

[cpp]  view plain copy
  1. View Code   
  2.   
  3. function out = adapthisteq(varargin)  
  4. %ADAPTHISTEQ Contrast-limited Adaptive Histogram Equalization (CLAHE).  
  5. %   ADAPTHISTEQ enhances the contrast of images by transforming the  
  6. %   values in the intensity image I.  Unlike HISTEQ, it operates on small  
  7. %   data regions (tiles), rather than the entire image. Each tile's   
  8. %   contrast is enhanced, so that the histogram of the output region  
  9. %   approximately matches the specified histogram. The neighboring tiles   
  10. %   are then combined using bilinear interpolation in order to eliminate  
  11. %   artificially induced boundaries.  The contrast, especially  
  12. %   in homogeneous areas, can be limited in order to avoid amplifying the  
  13. %   noise which might be present in the image.  
  14. %  
  15. %   J = ADAPTHISTEQ(I) Performs CLAHE on the intensity image I.  
  16. %  
  17. %   J = ADAPTHISTEQ(I,PARAM1,VAL1,PARAM2,VAL2...) sets various parameters.  
  18. %   Parameter names can be abbreviated, and case does not matter. Each   
  19. %   string parameter is followed by a value as indicated below:  
  20. %  
  21. %   'NumTiles'     Two-element vector of positive integers: [M N].  
  22. %                  [M N] specifies the number of tile rows and  
  23. %                  columns.  Both M and N must be at least 2.   
  24. %                  The total number of image tiles is equal to M*N.  
  25. %  
  26. %                  Default: [8 8].  
  27. %  
  28. %   'ClipLimit'    Real scalar from 0 to 1.  
  29. %                  'ClipLimit' limits contrast enhancement. Higher numbers   
  30. %                  result in more contrast.   
  31. %         
  32. %                  Default: 0.01.  
  33. %  
  34. %   'NBins'        Positive integer scalar.  
  35. %                  Sets number of bins for the histogram used in building a  
  36. %                  contrast enhancing transformation. Higher values result   
  37. %                  in greater dynamic range at the cost of slower processing  
  38. %                  speed.  
  39. %  
  40. %                  Default: 256.  
  41. %  
  42. %   'Range'        One of the strings: 'original' or 'full'.  
  43. %                  Controls the range of the output image data. If 'Range'   
  44. %                  is set to 'original', the range is limited to   
  45. %                  [min(I(:)) max(I(:))]. Otherwise, by default, or when   
  46. %                  'Range' is set to 'full', the full range of the output   
  47. %                  image class is used (e.g. [0 255] for uint8).  
  48. %  
  49. %                  Default: 'full'.  
  50. %  
  51. %   'Distribution' Distribution can be one of three strings: 'uniform',  
  52. %                  'rayleigh''exponential'.  
  53. %                  Sets desired histogram shape for the image tiles, by   
  54. %                  specifying a distribution type.  
  55. %  
  56. %                  Default: 'uniform'.  
  57. %  
  58. %   'Alpha'        Nonnegative real scalar.  
  59. %                  'Alpha' is a distribution parameter, which can be supplied   
  60. %                  when 'Dist' is set to either 'rayleigh' or 'exponential'.  
  61. %  
  62. %                  Default: 0.4.  
  63. %  
  64. %   Notes  
  65. %   -----  
  66. %   - 'NumTiles' specify the number of rectangular contextual regions (tiles)  
  67. %     into which the image is divided. The contrast transform function is  
  68. %     calculated for each of these regions individually. The optimal number of  
  69. %     tiles depends on the type of the input image, and it is best determined  
  70. %     through experimentation.  
  71. %  
  72. %   - The 'ClipLimit' is a contrast factor that prevents over-saturation of the  
  73. %     image specifically in homogeneous areas.  These areas are characterized  
  74. %     by a high peak in the histogram of the particular image tile due to many  
  75. %     pixels falling inside the same gray level range. Without the clip limit,  
  76. %     the adaptive histogram equalization technique could produce results that,  
  77. %     in some cases, are worse than the original image.  
  78. %  
  79. %   - ADAPTHISTEQ can use Uniform, Rayleigh, or Exponential distribution as  
  80. %     the basis for creating the contrast transform function. The distribution  
  81. %     that should be used depends on the type of the input image.  
  82. %     For example, underwater imagery appears to look more natural when the  
  83. %     Rayleigh distribution is used.  
  84. %  
  85. %   Class Support  
  86. %   -------------  
  87. %   Intensity image I can be uint8, uint16, int16, double, or single.  
  88. %   The output image J has the same class as I.  
  89. %  
  90. %   Example 1  
  91. %   ---------  
  92. %   Apply Contrast-Limited Adaptive Histogram Equalization to an   
  93. %   image and display the results.   
  94. %  
  95. %      I = imread('tire.tif');  
  96. %      A = adapthisteq(I,'clipLimit',0.02,'Distribution','rayleigh');  
  97. %      figure, imshow(I);  
  98. %      figure, imshow(A);  
  99. %  
  100. %   Example 2  
  101. %   ---------  
  102. %    
  103. %   Apply Contrast-Limited Adaptive Histogram Equalization to a color   
  104. %   photograph.  
  105. %  
  106. %      [X MAP] = imread('shadow.tif');  
  107. %      RGB = ind2rgb(X,MAP); % convert indexed image to truecolor format  
  108. %      cform2lab = makecform('srgb2lab');  
  109. %      LAB = applycform(RGB, cform2lab); %convert image to L*a*b color space  
  110. %      L = LAB(:,:,1)/100; % scale the values to range from 0 to 1  
  111. %      LAB(:,:,1) = adapthisteq(L,'NumTiles',[8 8],'ClipLimit',0.005)*100;  
  112. %      cform2srgb = makecform('lab2srgb');  
  113. %      J = applycform(LAB, cform2srgb); %convert back to RGB  
  114. %      figure, imshow(RGB); %display the results  
  115. %      figure, imshow(J);  
  116. %  
  117. %   See also HISTEQ.  
  118.   
  119. %   Copyright 1993-2010 The MathWorks, Inc.  
  120. %   $Revision: 1.1.6.12 $  $Date: 2011/08/09 17:48:54 $  
  121.   
  122. %   References:  
  123. %      Karel Zuiderveld, "Contrast Limited Adaptive Histogram Equalization",  
  124. %      Graphics Gems IV, p. 474-485, code: p. 479-484  
  125. %  
  126. %      Hanumant Singh, Woods Hole Oceanographic Institution, personal  
  127. %      communication  
  128.   
  129. %--------------------------- The algorithm ----------------------------------  
  130. %  
  131. %  1. Obtain all the inputs:   
  132. %    * image  
  133. %    * number of regions in row and column directions  
  134. %    * number of bins for the histograms used in building image transform  
  135. %      function (dynamic range)  
  136. %    * clip limit for contrast limiting (normalized from 0 to 1)  
  137. %    * other miscellaneous options  
  138. %  2. Pre-process the inputs:    
  139. %    * determine real clip limit from the normalized value  
  140. %    * if necessary, pad the image before splitting it into regions  
  141. %  3. Process each contextual region (tile) thus producing gray level mappings  
  142. %    * extract a single image region  
  143. %    * make a histogram for this region using the specified number of bins  
  144. %    * clip the histogram using clip limit  
  145. %    * create a mapping (transformation function) for this region  
  146. %  4. Interpolate gray level mappings in order to assemble final CLAHE image  
  147. %    * extract cluster of four neighboring mapping functions  
  148. %    * process image region partly overlapping each of the mapping tiles  
  149. %    * extract a single pixel, apply four mappings to that pixel, and   
  150. %      interpolate between the results to obtain the output pixel; repeat  
  151. %      over the entire image  
  152. %  
  153. %  See code for further details.  
  154. %  
  155. %-----------------------------------------------------------------------------  
  156.   
  157. [I, selectedRange, fullRange, numTiles, dimTile, clipLimit, numBins, ...  
  158.  noPadRect, distribution, alpha, int16ClassChange] = parseInputs(varargin{:});  
  159.   
  160. tileMappings = makeTileMappings(I, numTiles, dimTile, numBins, clipLimit, ...  
  161.                                 selectedRange, fullRange, distribution, alpha);  
  162.   
  163. %Synthesize the output image based on the individual tile mappings.   
  164. out = makeClaheImage(I, tileMappings, numTiles, selectedRange, numBins,...  
  165.                      dimTile);  
  166.   
  167. if int16ClassChange  
  168.   % Change uint16 back to int16 so output has same class as input.  
  169.   out = uint16toint16(out);  
  170. end  
  171.   
  172. if ~isempty(noPadRect) %do we need to remove padding?  
  173.   out = out(noPadRect.ulRow:noPadRect.lrRow, ...  
  174.             noPadRect.ulCol:noPadRect.lrCol);  
  175. end  
  176.   
  177. %-----------------------------------------------------------------------------  
  178.   
  179. function tileMappings = ...  
  180.     makeTileMappings(I, numTiles, dimTile, numBins, clipLimit,...  
  181.                      selectedRange, fullRange, distribution, alpha)  
  182.   
  183. numPixInTile = prod(dimTile);  
  184.   
  185. tileMappings = cell(numTiles);  
  186.   
  187. % extract and process each tile  
  188. imgCol = 1;  
  189. for col=1:numTiles(2),  
  190.   imgRow = 1;  
  191.   for row=1:numTiles(1),  
  192.       
  193.     tile = I(imgRow:imgRow+dimTile(1)-1,imgCol:imgCol+dimTile(2)-1);  
  194.   
  195.     % for speed, call MEX file directly thus avoiding costly  
  196.     % input parsing of imhist  
  197.     tileHist = imhistc(tile, numBins, 1, fullRange(2));   
  198.       
  199.     tileHist = clipHistogram(tileHist, clipLimit, numBins);  
  200.       
  201.     tileMapping = makeMapping(tileHist, selectedRange, fullRange, ...  
  202.                               numPixInTile, distribution, alpha);  
  203.       
  204.     % assemble individual tile mappings by storing them in a cell array;  
  205.     tileMappings{row,col} = tileMapping;  
  206.   
  207.     imgRow = imgRow + dimTile(1);      
  208.   end  
  209.   imgCol = imgCol + dimTile(2); % move to the next column of tiles  
  210. end  
  211.   
  212. %-----------------------------------------------------------------------------  
  213. % Calculate the equalized lookup table (mapping) based on cumulating the input   
  214. % histogram.  Note: lookup table is rescaled in the selectedRange [Min..Max].  
  215.   
  216. function mapping = makeMapping(imgHist, selectedRange, fullRange, ...  
  217.                                numPixInTile, distribution, alpha)  
  218.   
  219. histSum = cumsum(imgHist);  
  220. valSpread  = selectedRange(2) - selectedRange(1);  
  221.   
  222. switch distribution  
  223.  case 'uniform',  
  224.   scale =  valSpread/numPixInTile;  
  225.   mapping = min(selectedRange(1) + histSum*scale,...  
  226.                 selectedRange(2)); %limit to max  
  227.     
  228.  case 'rayleigh', % suitable for underwater imagery  
  229.   % pdf = (t./alpha^2).*exp(-t.^2/(2*alpha^2))*U(t)  
  230.   % cdf = 1-exp(-t.^2./(2*alpha^2))  
  231.   hconst = 2*alpha^2;  
  232.   vmax = 1 - exp(-1/hconst);  
  233.   val = vmax*(histSum/numPixInTile);  
  234.   val(val>=1) = 1-eps; % avoid log(0)  
  235.   temp = sqrt(-hconst*log(1-val));  
  236.   mapping = min(selectedRange(1)+temp*valSpread,...  
  237.                 selectedRange(2)); %limit to max  
  238.     
  239.  case 'exponential',  
  240.   % pdf = alpha*exp(-alpha*t)*U(t)  
  241.   % cdf = 1-exp(-alpha*t)  
  242.   vmax = 1 - exp(-alpha);  
  243.   val = (vmax*histSum/numPixInTile);  
  244.   val(val>=1) = 1-eps;  
  245.   temp = -1/alpha*log(1-val);  
  246.   mapping = min(selectedRange(1)+temp*valSpread,selectedRange(2));  
  247.     
  248.  otherwise,  
  249.   error(message('images:adapthisteq:distributionType')) %should never get here  
  250.     
  251. end  
  252.   
  253. %rescale the result to be between 0 and 1 for later use by the GRAYXFORM   
  254. %private mex function  
  255. mapping = mapping/fullRange(2);  
  256.   
  257. %-----------------------------------------------------------------------------  
  258. % This function clips the histogram according to the clipLimit and  
  259. % redistributes clipped pixels across bins below the clipLimit  
  260.   
  261. function imgHist = clipHistogram(imgHist, clipLimit, numBins)  
  262.   
  263. % total number of pixels overflowing clip limit in each bin  
  264. totalExcess = sum(max(imgHist - clipLimit,0));    
  265.   
  266. % clip the histogram and redistribute the excess pixels in each bin  
  267. avgBinIncr = floor(totalExcess/numBins);  
  268. upperLimit = clipLimit - avgBinIncr; % bins larger than this will be  
  269.                                      % set to clipLimit  
  270.   
  271. this loop should speed up the operation by putting multiple pixels  
  272. % into the "obvious" places first  
  273. for k=1:numBins  
  274.   if imgHist(k) > clipLimit  
  275.     imgHist(k) = clipLimit;  
  276.   else  
  277.     if imgHist(k) > upperLimit % high bin count  
  278.       totalExcess = totalExcess - (clipLimit - imgHist(k));  
  279.       imgHist(k) = clipLimit;  
  280.     else  
  281.       totalExcess = totalExcess - avgBinIncr;  
  282.       imgHist(k) = imgHist(k) + avgBinIncr;        
  283.     end  
  284.   end  
  285. end  
  286.   
  287. this loops redistributes the remaining pixels, one pixel at a time  
  288. k = 1;  
  289. while (totalExcess ~= 0)  
  290.   %keep increasing the step as fewer and fewer pixels remain for  
  291.   %the redistribution (spread them evenly)  
  292.   stepSize = max(floor(numBins/totalExcess),1);  
  293.   for m=k:stepSize:numBins  
  294.     if imgHist(m) < clipLimit  
  295.       imgHist(m) = imgHist(m)+1;  
  296.       totalExcess = totalExcess - 1; %reduce excess  
  297.       if totalExcess == 0  
  298.         break;  
  299.       end  
  300.     end  
  301.   end  
  302.     
  303.   k = k+1; %prevent from always placing the pixels in bin #1  
  304.   if k > numBins % start over if numBins was reached  
  305.     k = 1;  
  306.   end  
  307. end  
  308.   
  309. %-----------------------------------------------------------------------------  
  310. % This function interpolates between neighboring tile mappings to produce a   
  311. new mapping in order to remove artificially induced tile borders.    
  312. % Otherwise, these borders would become quite visible.  The resulting  
  313. % mapping is applied to the input image thus producing a CLAHE processed  
  314. % image.  
  315.   
  316. function claheI = makeClaheImage(I, tileMappings, numTiles, selectedRange,...  
  317.                                  numBins, dimTile)  
  318.   
  319. %initialize the output image to zeros (preserve the class of the input image)  
  320. claheI = I;  
  321. claheI(:) = 0;  
  322.   
  323. %compute the LUT for looking up original image values in the tile mappings,  
  324. %which we created earlier  
  325. if ~(isa(I,'double') || isa(I,'single'))  
  326.   k = selectedRange(1)+1 : selectedRange(2)+1;  
  327.   aLut = zeros(length(k),1);  
  328.   aLut(k) = (k-1)-selectedRange(1);  
  329.   aLut = aLut/(selectedRange(2)-selectedRange(1));  
  330. else  
  331.   % remap from 0..1 to 0..numBins-1  
  332.   if numBins ~= 1  
  333.     binStep = 1/(numBins-1);  
  334.     start = ceil(selectedRange(1)/binStep);  
  335.     stop  = floor(selectedRange(2)/binStep);  
  336.     k = start+1:stop+1;  
  337.     aLut(k) = 0:1/(length(k)-1):1;  
  338.   else  
  339.     aLut(1) = 0; %in case someone specifies numBins = 1, which is just silly  
  340.   end  
  341. end  
  342.   
  343. imgTileRow=1;  
  344. for k=1:numTiles(1)+1  
  345.   if k == 1  %special case: top row  
  346.     imgTileNumRows = dimTile(1)/2; %always divisible by 2 because of padding  
  347.     mapTileRows = [1 1];  
  348.   else   
  349.     if k == numTiles(1)+1 %special case: bottom row        
  350.       imgTileNumRows = dimTile(1)/2;  
  351.       mapTileRows = [numTiles(1) numTiles(1)];  
  352.     else %default values  
  353.       imgTileNumRows = dimTile(1);   
  354.       mapTileRows = [k-1, k]; %[upperRow lowerRow]  
  355.     end  
  356.   end  
  357.     
  358.   % loop over columns of the tileMappings cell array  
  359.   imgTileCol=1;  
  360.   for l=1:numTiles(2)+1  
  361.     if l == 1 %special case: left column  
  362.       imgTileNumCols = dimTile(2)/2;  
  363.       mapTileCols = [1, 1];  
  364.     else  
  365.       if l == numTiles(2)+1 % special case: right column  
  366.         imgTileNumCols = dimTile(2)/2;  
  367.         mapTileCols = [numTiles(2), numTiles(2)];  
  368.       else %default values  
  369.         imgTileNumCols = dimTile(2);  
  370.         mapTileCols = [l-1, l]; % right left  
  371.       end  
  372.     end  
  373.       
  374.     % Extract four tile mappings  
  375.     ulMapTile = tileMappings{mapTileRows(1), mapTileCols(1)};  
  376.     urMapTile = tileMappings{mapTileRows(1), mapTileCols(2)};  
  377.     blMapTile = tileMappings{mapTileRows(2), mapTileCols(1)};  
  378.     brMapTile = tileMappings{mapTileRows(2), mapTileCols(2)};  
  379.   
  380.     % Calculate the new greylevel assignments of pixels   
  381.     % within a submatrix of the image specified by imgTileIdx. This   
  382.     % is done by a bilinear interpolation between four different mappings   
  383.     % in order to eliminate boundary artifacts.  
  384.       
  385.     normFactor = imgTileNumRows*imgTileNumCols; %normalization factor    
  386.     imgTileIdx = {imgTileRow:imgTileRow+imgTileNumRows-1, ...  
  387.                  imgTileCol:imgTileCol+imgTileNumCols-1};  
  388.   
  389.     imgPixVals = grayxform(I(imgTileIdx{1},imgTileIdx{2}), aLut);  
  390.       
  391.     % calculate the weights used for linear interpolation between the  
  392.     % four mappings  
  393.     rowW = repmat((0:imgTileNumRows-1)',1,imgTileNumCols);  
  394.     colW = repmat(0:imgTileNumCols-1,imgTileNumRows,1);  
  395.     rowRevW = repmat((imgTileNumRows:-1:1)',1,imgTileNumCols);  
  396.     colRevW = repmat(imgTileNumCols:-1:1,imgTileNumRows,1);  
  397.       
  398.     claheI(imgTileIdx{1}, imgTileIdx{2}) = ...  
  399.         (rowRevW .* (colRevW .* double(grayxform(imgPixVals,ulMapTile)) + ...  
  400.                      colW    .* double(grayxform(imgPixVals,urMapTile)))+ ...  
  401.          rowW    .* (colRevW .* double(grayxform(imgPixVals,blMapTile)) + ...  
  402.                      colW    .* double(grayxform(imgPixVals,brMapTile))))...  
  403.         /normFactor;  
  404.       
  405.     imgTileCol = imgTileCol + imgTileNumCols;      
  406.   end %over tile cols  
  407.   imgTileRow = imgTileRow + imgTileNumRows;  
  408. end %over tile rows  
  409.   
  410. %-----------------------------------------------------------------------------  
  411.   
  412. function [I, selectedRange, fullRange, numTiles, dimTile, clipLimit,...  
  413.           numBins, noPadRect, distribution, alpha, ...  
  414.           int16ClassChange] = parseInputs(varargin)  
  415.   
  416. narginchk(1,13);  
  417.   
  418. I = varargin{1};  
  419. validateattributes(I, {'uint8''uint16''double''int16''single'}, ...  
  420.               {'real''2d''nonsparse''nonempty'}, ...  
  421.               mfilename, 'I', 1);  
  422.   
  423. % convert int16 to uint16  
  424. if isa(I,'int16')  
  425.   I = int16touint16(I);  
  426.   int16ClassChange = true;  
  427. else  
  428.   int16ClassChange = false;  
  429. end  
  430.   
  431. if any(size(I) < 2)  
  432.   error(message('images:adapthisteq:inputImageTooSmall'))  
  433. end  
  434.   
  435. %Other options  
  436. %%%%%%%%%%%%%%  
  437.   
  438. %Set the defaults  
  439. distribution = 'uniform';  
  440. alpha   = 0.4;  
  441.   
  442. if isa(I, 'double') || isa(I,'single')  
  443.   fullRange = [0 1];  
  444. else  
  445.   fullRange(1) = I(1);         %copy class of the input image  
  446.   fullRange(1:2) = [-Inf Inf]; %will be clipped to min and max  
  447.   fullRange = double(fullRange);  
  448. end  
  449.   
  450. selectedRange   = fullRange;  
  451.   
  452. %Set the default to 256 bins regardless of the data type;  
  453. %the user can override this value at any time  
  454. numBins = 256;  
  455. normClipLimit = 0.01;  
  456. numTiles = [8 8];  
  457.   
  458. checkAlpha = false;  
  459.   
  460. validStrings = {'NumTiles','ClipLimit','NBins','Distribution',...  
  461.                 'Alpha','Range'};  
  462.   
  463. if nargin > 1  
  464.   done = false;  
  465.   
  466.   idx = 2;  
  467.   while ~done  
  468.     input = varargin{idx};  
  469.     inputStr = validatestring(input, validStrings,mfilename,'PARAM',idx);  
  470.   
  471.     idx = idx+1; %advance index to point to the VAL portion of the input   
  472.   
  473.     if idx > nargin  
  474.       error(message('images:adapthisteq:missingValue', inputStr))  
  475.     end  
  476.       
  477.     switch inputStr  
  478.   
  479.      case 'NumTiles'  
  480.        numTiles = varargin{idx};  
  481.        validateattributes(numTiles, {'double'}, {'real''vector', ...  
  482.                            'integer''finite','positive'},...  
  483.                      mfilename, inputStr, idx);  
  484.   
  485.        if (any(size(numTiles) ~= [1,2]))  
  486.          error(message('images:adapthisteq:invalidNumTilesVector', inputStr))  
  487.        end  
  488.          
  489.        if any(numTiles < 2)  
  490.          error(message('images:adapthisteq:invalidNumTilesValue', inputStr))  
  491.        end  
  492.         
  493.      case 'ClipLimit'  
  494.       normClipLimit = varargin{idx};  
  495.       validateattributes(normClipLimit, {'double'}, ...  
  496.                     {'scalar','real','nonnegative'},...  
  497.                     mfilename, inputStr, idx);  
  498.         
  499.       if normClipLimit > 1  
  500.         error(message('images:adapthisteq:invalidClipLimit', inputStr))  
  501.       end  
  502.        
  503.      case 'NBins'  
  504.       numBins = varargin{idx};        
  505.       validateattributes(numBins, {'double'}, {'scalar','real','integer',...  
  506.                           'positive'}, mfilename, 'NBins', idx);  
  507.        
  508.      case 'Distribution'  
  509.       validDist = {'rayleigh','exponential','uniform'};  
  510.       distribution = validatestring(varargin{idx}, validDist, mfilename,...  
  511.                                   'Distribution', idx);  
  512.        
  513.      case 'Alpha'  
  514.       alpha = varargin{idx};  
  515.       validateattributes(alpha, {'double'},{'scalar','real',...  
  516.                           'nonnan','positive','finite'},...  
  517.                     mfilename, 'Alpha',idx);  
  518.       checkAlpha = true;  
  519.   
  520.      case 'Range'  
  521.       validRangeStrings = {'original','full'};  
  522.       rangeStr = validatestring(varargin{idx}, validRangeStrings,mfilename,...  
  523.                               'Range',idx);  
  524.         
  525.       if strmatch(rangeStr,'original')  
  526.         selectedRange = double([min(I(:)), max(I(:))]);  
  527.       end  
  528.        
  529.      otherwise  
  530.       error(message('images:adapthisteq:inputString')) %should never get here  
  531.     end  
  532.       
  533.     if idx >= nargin  
  534.       done = true;  
  535.     end  
  536.       
  537.     idx=idx+1;  
  538.   end  
  539. end  
  540.   
  541.   
  542. %% Pre-process the inputs  
  543. %%%%%%%%%%%%%%%%%%%%%%%%%%  
  544.   
  545. dimI = size(I);  
  546. dimTile = dimI ./ numTiles;  
  547.   
  548. %check if tile size is reasonable  
  549. if any(dimTile < 1)  
  550.   error(message('images:adapthisteq:inputImageTooSmallToSplit', num2str( numTiles )))  
  551. end  
  552.   
  553. if checkAlpha  
  554.   if strcmp(distribution,'uniform')  
  555.     error(message('images:adapthisteq:alphaShouldNotBeSpecified', distribution))  
  556.   end  
  557. end  
  558.   
  559. %check if the image needs to be padded; pad if necessary;  
  560. %padding occurs if any dimension of a single tile is an odd number  
  561. %and/or when image dimensions are not divisible by the selected   
  562. %number of tiles  
  563. rowDiv  = mod(dimI(1),numTiles(1)) == 0;  
  564. colDiv  = mod(dimI(2),numTiles(2)) == 0;  
  565.   
  566. if rowDiv && colDiv  
  567.   rowEven = mod(dimTile(1),2) == 0;  
  568.   colEven = mod(dimTile(2),2) == 0;    
  569. end  
  570.   
  571. noPadRect = [];  
  572. if  ~(rowDiv && colDiv && rowEven && colEven)  
  573.   padRow = 0;  
  574.   padCol = 0;  
  575.     
  576.   if ~rowDiv  
  577.     rowTileDim = floor(dimI(1)/numTiles(1)) + 1;  
  578.     padRow = rowTileDim*numTiles(1) - dimI(1);  
  579.   else  
  580.     rowTileDim = dimI(1)/numTiles(1);  
  581.   end  
  582.     
  583.   if ~colDiv  
  584.     colTileDim = floor(dimI(2)/numTiles(2)) + 1;  
  585.     padCol = colTileDim*numTiles(2) - dimI(2);  
  586.   else  
  587.     colTileDim = dimI(2)/numTiles(2);  
  588.   end  
  589.     
  590.   %check if tile dimensions are even numbers  
  591.   rowEven = mod(rowTileDim,2) == 0;  
  592.   colEven = mod(colTileDim,2) == 0;  
  593.     
  594.   if ~rowEven  
  595.     padRow = padRow+numTiles(1);  
  596.   end  
  597.   
  598.   if ~colEven  
  599.     padCol = padCol+numTiles(2);  
  600.   end  
  601.     
  602.   padRowPre  = floor(padRow/2);  
  603.   padRowPost = ceil(padRow/2);  
  604.   padColPre  = floor(padCol/2);  
  605.   padColPost = ceil(padCol/2);  
  606.     
  607.   I = padarray(I,[padRowPre  padColPre ],'symmetric','pre');  
  608.   I = padarray(I,[padRowPost padColPost],'symmetric','post');  
  609.   
  610.   %UL corner (Row, Col), LR corner (Row, Col)  
  611.   noPadRect.ulRow = padRowPre+1;  
  612.   noPadRect.ulCol = padColPre+1;  
  613.   noPadRect.lrRow = padRowPre+dimI(1);  
  614.   noPadRect.lrCol = padColPre+dimI(2);  
  615. end  
  616.   
  617. %redefine this variable to include the padding  
  618. dimI = size(I);  
  619.   
  620. %size of the single tile  
  621. dimTile = dimI ./ numTiles;  
  622.   
  623. %compute actual clip limit from the normalized value entered by the user  
  624. %maximum value of normClipLimit=1 results in standard AHE, i.e. no clipping;  
  625. %the minimum value minClipLimit would uniformly distribute the image pixels  
  626. %across the entire histogram, which would result in the lowest possible  
  627. %contrast value  
  628. numPixInTile = prod(dimTile);  
  629. minClipLimit = ceil(numPixInTile/numBins);  
  630. clipLimit = minClipLimit + round(normClipLimit*(numPixInTile-minClipLimit));  
  631.   
  632. %-----------------------------------------------------------------------------  


      参考上述代码,我分别用VB和C#实现了该算法,提供个编译好的文件给有兴趣研究该算法的朋友看看效果(不提供源代码的):

      http://files.cnblogs.com/Imageshop/CLAHE.rar

      截面如下:

   其中AHE算法可以认为是裁剪限幅为1的CLAHE算法,CLHE是水平网格和垂直网格都为1的算法。

   均衡分布方式和ALPHA的解释可参考matlab的代码.

   CLAHE算法很多时候比直接的直方图均衡化算法的效果要好很多,比如下面这幅图像:

      

  普通的直方图均衡化算法的效果为:

 

  而CALHE的效果如下:

    可以看出,在图像的上部,CALHE有效的抑制了噪音的增强。

    对于彩色图像,matlab的那个函数不支持,我这里采用了两种方式处理,一种是各通道分开处理,一种是每个通道都放在一起,对于彩色的图像似乎放在一起处理的效果要比分开处理好很多。 

     比如界面中那个图像,如果各通道放在一起处理,效果如下:

    

    而通道分开处理,则各通道的颜色不很匹配:

    

     对于处理速度,这个函数由于只有一些分开+插值计算,速度很快。如果图像不能倍分成整块,一般需要扩充边缘,matlab的处理方式是上下左右四条边对称镜像扩充。这种方式很合理。

     实例工程是用VB6编写的,由于VB不支持指针,在速度比C#之类的语言大概慢了50%左右。但是也还是很快的了。


  • 1
    点赞
  • 0
    评论
  • 4
    收藏
  • 一键三连
    一键三连
  • 扫一扫,分享海报

©️2021 CSDN 皮肤主题: 大白 设计师:CSDN官方博客 返回首页
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、C币套餐、付费专栏及课程。

余额充值