Self-adaptive Histogram equalization (AHE) and constrained contrast (CLAHE)

This article is translated from https://en.wikipedia.org/wiki/Adaptive_histogram_equalization, if there are errors also hope Haihan.

Adaptive histogram Equalization (Adaptive histogram equalization)

1. Introduction to Algorithms

Ahe is an image processing technology to improve image contrast, compared with traditional (normal) histogram equalization, the difference is mainly that ahe by calculating the histogram of each significant area of the image to redistribute the luminance value of the image, so it is more suitable for improving the local contrast of the image, and enhance image edge information, facilitate segmentation.

However, Ahe has a flaw that he enhances contrast while also enhancing the noise of the homogeneous (uniform) area of the image, so as an improvement of ahe, Clahe can effectively reduce the noise enhancement.

2. The meaning of the algorithm

The traditional histogram equalization uses the same histogram transform when the whole image is transformed, this method is very good for those images with uniform pixel distribution, but it can not play a significant effect on the image which contains the obviously lighter or darker area.

Ahe performs histogram enhancement on each pixel by calculating the transformation function of each pixel neighborhood, which was first used in the display image enhancement of the spacecraft cockpit. The simplest form is to equalize each pixel based on the histogram of the square neighborhood of the pixel, such as the histogram-based enhancement, which is exactly the same as the normal histogram enhancement, since the transformation function is proportionate to the cumulative distribution function (CDF) of the Pixel neighborhood.

The boundary pixels of the image should be handled specially because the neighbors of these pixels are not fully contained within the image, such as the blue box, which can be resolved by mirroring the pixel rows or columns, and it is undesirable to replicate directly to the boundary pixel rows, as this results in a high spike in the pixel neighborhood histogram.

The nature of the 3.AHE algorithm

(1) The size of the neighborhood is a scale parameter of the neighborhood length, when the parameter is large, the contrast is reduced; Conversely, the contrast is enhanced;

(2) Due to the nature of histogram enhancement, the result of pixel transformation in Ahe is proportional to its order in the neighborhood pixels, so we can use specialized hardware to realize the comparison between this neighborhood center pixel and other pixels, and a non-normalized method of calculation, This can be achieved by adding a value smaller than the center pixel to each pixel, or by adding a neighborhood pixel mean on each pixel.

(3) If a pixel neighborhood is homogeneous or homogeneous, then the histogram of the neighborhood may be sharp, that is, a very high peak, then the transformation function will be a very narrow range of pixel values mapped to the full range of pixel values of the transformation results, which is why the above mentioned ahe to enhance the uniform pixel neighborhood noise reasons.

Adaptive histogram enhancement with limited contrast (contrast Limited AHE)

1. Fundamentals

The biggest difference between Clahe and Ahe is that the former limits the contrast, which can also be applied to the global histogram equalization, that is, contrast Limited he, abbreviated to CLHE, but it is seldom used in practice. In Clahe, each pixel neighborhood is limited in contrast to obtain the corresponding transformation function, which is used to reduce the noise enhancement in the Ahe, which is mainly achieved by limiting the contrast enhancement in the ahe. The enhancement of neighborhood noise around pixels is mainly caused by the slope of the transformation function, because the noise of the pixel neighborhood is proportional to the neighborhood CDF, so it is also proportional to the value of the neighborhood histogram in the center pixel position, and the Clahe can limit the contrast. Because it trims the histogram at the specified threshold before computing the CDF of the neighborhood, as shown, this approach not only limits the slope of the CDF, but also limits the slope of the transformation function, where the threshold used to cut the histogram is called the clipping limit system (clip limit), This parameter depends not only on the normalization of the histogram, but also on the size of the pixel neighborhood, which is usually set from 3 to 4.

, it does not throw away the trimmed portion of the histogram directly, but instead uses a more appropriate approach, which is to redistribute the trimmed portions into the histogram to generate a new histogram.

Also, in this redistribution, it is possible that some areas will again exceed clip limit, such as the area labeled Green in the image on the right, when we need to adjust the clip limit to achieve the actual need, if still can not let you satisfied, You can try the recursive procedure to repeat the step of redistribution of the histogram until the portion of the excess affects our results negligible.

2. Interpolation acceleration

The Adaptive histogram equalization algorithm mentioned above, whether or not the contrast is limited, requires that the neighborhood histogram and corresponding transformation function be computed for the pixel neighborhood of each pixel, a process that is quite time-consuming the.

Therefore, on the basis of guaranteeing the effect of the algorithm, the interpolation algorithm is used to improve the computational efficiency, as follows:

-A rectangular tile that divides the image into the same size, such as right, is usually divided into 8x8, i.e. 64 tiles;

-Each tile's histogram, CDF, and transform functions are computed separately, and for the center pixel of the tiles (the black rectangular block on the left), it is possible to directly use the transform function of the tile, because it satisfies the transformation function of the above calculation. Other pixels need to be interpolated by a transform function with at least four pixels nearest to it, which is the tiles (i.e., the neighborhood domain, which can also be 8 neighbors). Where the pixels in the left Middle blue area are marked with bilinear interpolation (note that this part of the pixel occupies a large portion of the image); The pixel points of the boundary part (the area of the green mark) are worth the linear interpolation, while the pixels near the corner of the image (red mark) are obtained by the transform function of the tile at the corner point.

Where the interpolation coefficients reflect the position between the center pixels of the nearest tiles, the interpolation process is continuous because the interpolated pixel points can be approximated as a central pixel.

Although this process increases the computational amount of linear interpolation, it greatly reduces the number of calculation of the transformation function, so it improves the computational speed.

Source code reference for 3.CLAHE

/* * ANSI C code from the article * "Contrast Limited Adaptive histogram equalization" * by Karel Zuiderveld, [email  Protected] * in "Graphics Gems IV", Academic Press, 1994 * * * These functions implement contrast Limited Adaptive Histog Ram equalization.  * The main routine (CLAHE) expects an input image which is stored contiguously in * memory; The CLAHE output image overwrites the original input image and have the * same minimum and maximum values (which must be P Rovided by the user).  * This implementation assumes that the X-and Y image resolutions is 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-contextual regions is not recommended. * * The code is ANSI-C and is AlSo C + + compliant.  * * Author:karel Zuiderveld, computer Vision Group, * Utrecht, the Netherlands ([email protected])     */#ifdef byte_imagetypedef unsigned char kz_pixel_t;     /* for 8 bit-per-pixel images */#define Uinr_of_grey (#elsetypedef) unsigned short kz_pixel_t; /* for 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, Unsigne d int UINRX, unsigned int uinry, unsigned int uinrbins, float fcliplimit);/*********************** 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*); St          atic 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 Interpol Ate (kz_pixel_t*, int, unsigned long*, unsigned long*, unsigned long*, unsigned long*, unsigned int, unsigned int, kz_p ixel_t*)/************** Start of actual code **************/#include <stdlib.h>/* to get prototype      S 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_pi xel_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-i Mage 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 Max Imum of output image) * Uinrx-number of contextial regions in the X direction (min 2, max uimax_reg_x) * Uinry-num ber of contextial regions in the Y direction (min 2, max uimax_reg_y) * Uinrbins-number of Greybins for histogram ("dy Namic range ") * Float fcliplimit-normalized cliplimit (higher values give more contrast) * The number of" effective "G Reylevels in the output image was set by Uinrbins; selecting * A small value (eg.) speeds up processing and still produce an output image of * good quality. The output image would 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;       /* # Regions X-direction too large */if (Uinry > Uimax_reg_y) return-2;      /* # 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 least4 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 + = u        ixsize) {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 >& Gt  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 is 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;    Long lbinexcess;  ulnrexcess = 0;    Pulbinpointer = Pulhistogram; for (i = 0; i < uinrgreylevels; i++) {/* Calculate total number of excess pixels */lbinexcess = (long) pulbinpoint    Er[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-= Pul Histogram[i]-ulupper;      Pulhistogram[i]=ulcliplimit;      } else {/* Low bin count */ulnrexcess-= ULBININCR; Pulhistogram[i] + = ULBININCR; }}} while (ulnrexcess) {/* Redistribute remaining excess */Pulendpointer = &pulhistogram[uinrgr Eylevels];    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 */}}}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 histogr Am. 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 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, u nsigned long Ulnrofpixels)/* This function calculates the equalized lookup table (mapping) by * cumulating the input Histo Gram. 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-p Ointer to Input/output image * Uixres-resolution of image in X-direction * pulmap*-mappings of Greylevels fro M histograms * Uixsize-uixsize of Image Submatrix * Uiysize-uiysize of Image Submatrix * Plut-look  Up 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 do 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; u nsigned 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 both, use division */for (Uiycoef = 0, Uiyinvcoef = Uiysiz E     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); }    }    }}

Self-adaptive Histogram equalization (AHE) and constrained contrast (CLAHE)