/* * ANSI C code from the article * "Contrast Limited Adaptive Histogram Equalization" * by Karel Zuiderveld, karel@cv.ruu.nl * in "Graphics Gems IV", Academic Press, 1994 * * *  These functions implement Contrast Limited Adaptive Histogram Equalization. *  The main routine (CLAHE) expects an input image that is stored contiguously in *  memory;  the CLAHE output image overwrites the original input image and has the *  same minimum and maximum values (which must be provided by the user). *  This implementation assumes that the X- and Y image resolutions are an integer *  multiple of the X- and Y sizes of the contextual regions. A check on various other *  error conditions is performed. * *  #define the symbol BYTE_IMAGE to make this implementation suitable for *  8-bit images. The maximum number of contextual regions can be redefined *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256 *  contextual regions is not recommended. * *  The code is ANSI-C and is also C++ compliant. * *  Author: Karel Zuiderveld, Computer Vision Research Group, *         Utrecht, The Netherlands (karel@cv.ruu.nl) */#ifdef BYTE_IMAGEtypedef unsigned char kz_pixel_t;     /* for 8 bit-per-pixel images */#define uiNR_OF_GREY (256)#elsetypedef unsigned short kz_pixel_t;     /* for 12 bit-per-pixel images (default) */# define uiNR_OF_GREY (4096)#endif/******** Prototype of CLAHE function. Put this in a separate include file. *****/int CLAHE(kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes, kz_pixel_t Min,      kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,      unsigned int uiNrBins, float fCliplimit);/*********************** Local prototypes ************************/static void ClipHistogram (unsigned long*, unsigned int, unsigned long);static void MakeHistogram (kz_pixel_t*, unsigned int, unsigned int, unsigned int,        unsigned long*, unsigned int, kz_pixel_t*);static void MapHistogram (unsigned long*, kz_pixel_t, kz_pixel_t,           unsigned int, unsigned long);static void MakeLut (kz_pixel_t*, kz_pixel_t, kz_pixel_t, unsigned int);static void Interpolate (kz_pixel_t*, int, unsigned long*, unsigned long*,    unsigned long*, unsigned long*, unsigned int, unsigned int, kz_pixel_t*);/**************     Start of actual code **************/#include <stdlib.h>             /* To get prototypes of malloc() and free() */const unsigned int uiMAX_REG_X = 16;      /* max. # contextual regions in x-direction */const unsigned int uiMAX_REG_Y = 16;      /* max. # contextual regions in y-direction *//************************** main function CLAHE ******************/int CLAHE (kz_pixel_t* pImage, unsigned int uiXRes, unsigned int uiYRes,     kz_pixel_t Min, kz_pixel_t Max, unsigned int uiNrX, unsigned int uiNrY,          unsigned int uiNrBins, float fCliplimit)/*   pImage - Pointer to the input/output image *   uiXRes - Image resolution in the X direction *   uiYRes - Image resolution in the Y direction *   Min - Minimum greyvalue of input image (also becomes minimum of output image) *   Max - Maximum greyvalue of input image (also becomes maximum of output image) *   uiNrX - Number of contextial regions in the X direction (min 2, max uiMAX_REG_X) *   uiNrY - Number of contextial regions in the Y direction (min 2, max uiMAX_REG_Y) *   uiNrBins - Number of greybins for histogram ("dynamic range") *   float fCliplimit - Normalized cliplimit (higher values give more contrast) * The number of "effective" greylevels in the output image is set by uiNrBins; selecting * a small value (eg. 128) speeds up processing and still produce an output image of * good quality. The output image will have the same minimum and maximum value as the input * image. A clip limit smaller than 1 results in standard (non-contrast limited) AHE. */{    unsigned int uiX, uiY;          /* counters */    unsigned int uiXSize, uiYSize, uiSubX, uiSubY; /* size of context. reg. and subimages */    unsigned int uiXL, uiXR, uiYU, uiYB;  /* auxiliary variables interpolation routine */    unsigned long ulClipLimit, ulNrPixels;/* clip limit and region pixel count */    kz_pixel_t* pImPointer;           /* pointer to image */    kz_pixel_t aLUT[uiNR_OF_GREY];        /* lookup table used for scaling of input image */    unsigned long* pulHist, *pulMapArray; /* pointer to histogram and mappings*/    unsigned long* pulLU, *pulLB, *pulRU, *pulRB; /* auxiliary pointers interpolation */    if (uiNrX > uiMAX_REG_X) return -1;       /* # of regions x-direction too large */    if (uiNrY > uiMAX_REG_Y) return -2;       /* # of regions y-direction too large */    if (uiXRes % uiNrX) return -3;      /* x-resolution no multiple of uiNrX */    if (uiYRes & uiNrY) return -4;      /* y-resolution no multiple of uiNrY */    if (Max >= uiNR_OF_GREY) return -5;       /* maximum too large */    if (Min >= Max) return -6;          /* minimum equal or larger than maximum */    if (uiNrX < 2 || uiNrY < 2) return -7;/* at least 4 contextual regions required */    if (fCliplimit == 1.0) return 0;      /* is OK, immediately returns original image. */    if (uiNrBins == 0) uiNrBins = 128;      /* default value when not specified */    pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);    if (pulMapArray == 0) return -8;      /* Not enough memory! (try reducing uiNrBins) */    uiXSize = uiXRes/uiNrX; uiYSize = uiYRes/uiNrY;  /* Actual size of contextual regions */    ulNrPixels = (unsigned long)uiXSize * (unsigned long)uiYSize;    if(fCliplimit > 0.0) {          /* Calculate actual cliplimit     */       ulClipLimit = (unsigned long) (fCliplimit * (uiXSize * uiYSize) / uiNrBins);       ulClipLimit = (ulClipLimit < 1UL) ? 1UL : ulClipLimit;    }    else ulClipLimit = 1UL<<14;          /* Large value, do not clip (AHE) */    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;    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 */    }    }}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);        }    }    }}

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

  1 function out = adapthisteq(varargin)  2 %ADAPTHISTEQ Contrast-limited Adaptive Histogram Equalization (CLAHE).  3 %   ADAPTHISTEQ enhances the contrast of images by transforming the  4 %   values in the intensity image I.  Unlike HISTEQ, it operates on small  5 %   data regions (tiles), rather than the entire image. Each tile‘s   6 %   contrast is enhanced, so that the histogram of the output region  7 %   approximately matches the specified histogram. The neighboring tiles   8 %   are then combined using bilinear interpolation in order to eliminate  9 %   artificially induced boundaries.  The contrast, especially 10 %   in homogeneous areas, can be limited in order to avoid amplifying the 11 %   noise which might be present in the image. 12 % 13 %   J = ADAPTHISTEQ(I) Performs CLAHE on the intensity image I. 14 % 15 %   J = ADAPTHISTEQ(I,PARAM1,VAL1,PARAM2,VAL2...) sets various parameters. 16 %   Parameter names can be abbreviated, and case does not matter. Each  17 %   string parameter is followed by a value as indicated below: 18 % 19 %   ‘NumTiles‘     Two-element vector of positive integers: [M N]. 20 %                  [M N] specifies the number of tile rows and 21 %                  columns.  Both M and N must be at least 2.  22 %                  The total number of image tiles is equal to M*N. 23 % 24 %                  Default: [8 8]. 25 % 26 %   ‘ClipLimit‘    Real scalar from 0 to 1. 27 %                  ‘ClipLimit‘ limits contrast enhancement. Higher numbers  28 %                  result in more contrast.  29 %        30 %                  Default: 0.01. 31 % 32 %   ‘NBins‘        Positive integer scalar. 33 %                  Sets number of bins for the histogram used in building a 34 %                  contrast enhancing transformation. Higher values result  35 %                  in greater dynamic range at the cost of slower processing 36 %                  speed. 37 % 38 %                  Default: 256. 39 % 40 %   ‘Range‘        One of the strings: ‘original‘ or ‘full‘. 41 %                  Controls the range of the output image data. If ‘Range‘  42 %                  is set to ‘original‘, the range is limited to  43 %                  [min(I(:)) max(I(:))]. Otherwise, by default, or when  44 %                  ‘Range‘ is set to ‘full‘, the full range of the output  45 %                  image class is used (e.g. [0 255] for uint8). 46 % 47 %                  Default: ‘full‘. 48 % 49 %   ‘Distribution‘ Distribution can be one of three strings: ‘uniform‘, 50 %                  ‘rayleigh‘, ‘exponential‘. 51 %                  Sets desired histogram shape for the image tiles, by  52 %                  specifying a distribution type. 53 % 54 %                  Default: ‘uniform‘. 55 % 56 %   ‘Alpha‘        Nonnegative real scalar. 57 %                  ‘Alpha‘ is a distribution parameter, which can be supplied  58 %                  when ‘Dist‘ is set to either ‘rayleigh‘ or ‘exponential‘. 59 % 60 %                  Default: 0.4. 61 % 62 %   Notes 63 %   ----- 64 %   - ‘NumTiles‘ specify the number of rectangular contextual regions (tiles) 65 %     into which the image is divided. The contrast transform function is 66 %     calculated for each of these regions individually. The optimal number of 67 %     tiles depends on the type of the input image, and it is best determined 68 %     through experimentation. 69 % 70 %   - The ‘ClipLimit‘ is a contrast factor that prevents over-saturation of the 71 %     image specifically in homogeneous areas.  These areas are characterized 72 %     by a high peak in the histogram of the particular image tile due to many 73 %     pixels falling inside the same gray level range. Without the clip limit, 74 %     the adaptive histogram equalization technique could produce results that, 75 %     in some cases, are worse than the original image. 76 % 77 %   - ADAPTHISTEQ can use Uniform, Rayleigh, or Exponential distribution as 78 %     the basis for creating the contrast transform function. The distribution 79 %     that should be used depends on the type of the input image. 80 %     For example, underwater imagery appears to look more natural when the 81 %     Rayleigh distribution is used. 82 % 83 %   Class Support 84 %   ------------- 85 %   Intensity image I can be uint8, uint16, int16, double, or single. 86 %   The output image J has the same class as I. 87 % 88 %   Example 1 89 %   --------- 90 %   Apply Contrast-Limited Adaptive Histogram Equalization to an  91 %   image and display the results.  92 % 93 %      I = imread(‘tire.tif‘); 94 %      A = adapthisteq(I,‘clipLimit‘,0.02,‘Distribution‘,‘rayleigh‘); 95 %      figure, imshow(I); 96 %      figure, imshow(A); 97 % 98 %   Example 2 99 %   ---------100 %  101 %   Apply Contrast-Limited Adaptive Histogram Equalization to a color 102 %   photograph.103 %104 %      [X MAP] = imread(‘shadow.tif‘);105 %      RGB = ind2rgb(X,MAP); % convert indexed image to truecolor format106 %      cform2lab = makecform(‘srgb2lab‘);107 %      LAB = applycform(RGB, cform2lab); %convert image to L*a*b color space108 %      L = LAB(:,:,1)/100; % scale the values to range from 0 to 1109 %      LAB(:,:,1) = adapthisteq(L,‘NumTiles‘,[8 8],‘ClipLimit‘,0.005)*100;110 %      cform2srgb = makecform(‘lab2srgb‘);111 %      J = applycform(LAB, cform2srgb); %convert back to RGB112 %      figure, imshow(RGB); %display the results113 %      figure, imshow(J);114 %115 %   See also HISTEQ.116 117 %   Copyright 1993-2010 The MathWorks, Inc.118 %   $Revision: 1.1.6.12 $  $Date: 2011/08/09 17:48:54 $119 120 %   References:121 %      Karel Zuiderveld, "Contrast Limited Adaptive Histogram Equalization",122 %      Graphics Gems IV, p. 474-485, code: p. 479-484123 %124 %      Hanumant Singh, Woods Hole Oceanographic Institution, personal125 %      communication126 127 %--------------------------- The algorithm ----------------------------------128 %129 %  1. Obtain all the inputs: 130 %    * image131 %    * number of regions in row and column directions132 %    * number of bins for the histograms used in building image transform133 %      function (dynamic range)134 %    * clip limit for contrast limiting (normalized from 0 to 1)135 %    * other miscellaneous options136 %  2. Pre-process the inputs:  137 %    * determine real clip limit from the normalized value138 %    * if necessary, pad the image before splitting it into regions139 %  3. Process each contextual region (tile) thus producing gray level mappings140 %    * extract a single image region141 %    * make a histogram for this region using the specified number of bins142 %    * clip the histogram using clip limit143 %    * create a mapping (transformation function) for this region144 %  4. Interpolate gray level mappings in order to assemble final CLAHE image145 %    * extract cluster of four neighboring mapping functions146 %    * process image region partly overlapping each of the mapping tiles147 %    * extract a single pixel, apply four mappings to that pixel, and 148 %      interpolate between the results to obtain the output pixel; repeat149 %      over the entire image150 %151 %  See code for further details.152 %153 %-----------------------------------------------------------------------------154 155 [I, selectedRange, fullRange, numTiles, dimTile, clipLimit, numBins, ...156  noPadRect, distribution, alpha, int16ClassChange] = parseInputs(varargin{:});157 158 tileMappings = makeTileMappings(I, numTiles, dimTile, numBins, clipLimit, ...159                                 selectedRange, fullRange, distribution, alpha);160 161 %Synthesize the output image based on the individual tile mappings. 162 out = makeClaheImage(I, tileMappings, numTiles, selectedRange, numBins,...163                      dimTile);164 165 if int16ClassChange166   % Change uint16 back to int16 so output has same class as input.167   out = uint16toint16(out);168 end169 170 if ~isempty(noPadRect) %do we need to remove padding?171   out = out(noPadRect.ulRow:noPadRect.lrRow, ...172             noPadRect.ulCol:noPadRect.lrCol);173 end174 175 %-----------------------------------------------------------------------------176 177 function tileMappings = ...178     makeTileMappings(I, numTiles, dimTile, numBins, clipLimit,...179                      selectedRange, fullRange, distribution, alpha)180 181 numPixInTile = prod(dimTile);182 183 tileMappings = cell(numTiles);184 185 % extract and process each tile186 imgCol = 1;187 for col=1:numTiles(2),188   imgRow = 1;189   for row=1:numTiles(1),190     191     tile = I(imgRow:imgRow+dimTile(1)-1,imgCol:imgCol+dimTile(2)-1);192 193     % for speed, call MEX file directly thus avoiding costly194     % input parsing of imhist195     tileHist = imhistc(tile, numBins, 1, fullRange(2)); 196     197     tileHist = clipHistogram(tileHist, clipLimit, numBins);198     199     tileMapping = makeMapping(tileHist, selectedRange, fullRange, ...200                               numPixInTile, distribution, alpha);201     202     % assemble individual tile mappings by storing them in a cell array;203     tileMappings{row,col} = tileMapping;204 205     imgRow = imgRow + dimTile(1);    206   end207   imgCol = imgCol + dimTile(2); % move to the next column of tiles208 end209 210 %-----------------------------------------------------------------------------211 % Calculate the equalized lookup table (mapping) based on cumulating the input 212 % histogram.  Note: lookup table is rescaled in the selectedRange [Min..Max].213 214 function mapping = makeMapping(imgHist, selectedRange, fullRange, ...215                                numPixInTile, distribution, alpha)216 217 histSum = cumsum(imgHist);218 valSpread  = selectedRange(2) - selectedRange(1);219 220 switch distribution221  case ‘uniform‘,222   scale =  valSpread/numPixInTile;223   mapping = min(selectedRange(1) + histSum*scale,...224                 selectedRange(2)); %limit to max225   226  case ‘rayleigh‘, % suitable for underwater imagery227   % pdf = (t./alpha^2).*exp(-t.^2/(2*alpha^2))*U(t)228   % cdf = 1-exp(-t.^2./(2*alpha^2))229   hconst = 2*alpha^2;230   vmax = 1 - exp(-1/hconst);231   val = vmax*(histSum/numPixInTile);232   val(val>=1) = 1-eps; % avoid log(0)233   temp = sqrt(-hconst*log(1-val));234   mapping = min(selectedRange(1)+temp*valSpread,...235                 selectedRange(2)); %limit to max236   237  case ‘exponential‘,238   % pdf = alpha*exp(-alpha*t)*U(t)239   % cdf = 1-exp(-alpha*t)240   vmax = 1 - exp(-alpha);241   val = (vmax*histSum/numPixInTile);242   val(val>=1) = 1-eps;243   temp = -1/alpha*log(1-val);244   mapping = min(selectedRange(1)+temp*valSpread,selectedRange(2));245   246  otherwise,247   error(message(‘images:adapthisteq:distributionType‘)) %should never get here248   249 end250 251 %rescale the result to be between 0 and 1 for later use by the GRAYXFORM 252 %private mex function253 mapping = mapping/fullRange(2);254 255 %-----------------------------------------------------------------------------256 % This function clips the histogram according to the clipLimit and257 % redistributes clipped pixels across bins below the clipLimit258 259 function imgHist = clipHistogram(imgHist, clipLimit, numBins)260 261 % total number of pixels overflowing clip limit in each bin262 totalExcess = sum(max(imgHist - clipLimit,0));  263 264 % clip the histogram and redistribute the excess pixels in each bin265 avgBinIncr = floor(totalExcess/numBins);266 upperLimit = clipLimit - avgBinIncr; % bins larger than this will be267                                      % set to clipLimit268 269 % this loop should speed up the operation by putting multiple pixels270 % into the "obvious" places first271 for k=1:numBins272   if imgHist(k) > clipLimit273     imgHist(k) = clipLimit;274   else275     if imgHist(k) > upperLimit % high bin count276       totalExcess = totalExcess - (clipLimit - imgHist(k));277       imgHist(k) = clipLimit;278     else279       totalExcess = totalExcess - avgBinIncr;280       imgHist(k) = imgHist(k) + avgBinIncr;      281     end282   end283 end284 285 % this loops redistributes the remaining pixels, one pixel at a time286 k = 1;287 while (totalExcess ~= 0)288   %keep increasing the step as fewer and fewer pixels remain for289   %the redistribution (spread them evenly)290   stepSize = max(floor(numBins/totalExcess),1);291   for m=k:stepSize:numBins292     if imgHist(m) < clipLimit293       imgHist(m) = imgHist(m)+1;294       totalExcess = totalExcess - 1; %reduce excess295       if totalExcess == 0296         break;297       end298     end299   end300   301   k = k+1; %prevent from always placing the pixels in bin #1302   if k > numBins % start over if numBins was reached303     k = 1;304   end305 end306 307 %-----------------------------------------------------------------------------308 % This function interpolates between neighboring tile mappings to produce a 309 % new mapping in order to remove artificially induced tile borders.  310 % Otherwise, these borders would become quite visible.  The resulting311 % mapping is applied to the input image thus producing a CLAHE processed312 % image.313 314 function claheI = makeClaheImage(I, tileMappings, numTiles, selectedRange,...315                                  numBins, dimTile)316 317 %initialize the output image to zeros (preserve the class of the input image)318 claheI = I;319 claheI(:) = 0;320 321 %compute the LUT for looking up original image values in the tile mappings,322 %which we created earlier323 if ~(isa(I,‘double‘) || isa(I,‘single‘))324   k = selectedRange(1)+1 : selectedRange(2)+1;325   aLut = zeros(length(k),1);326   aLut(k) = (k-1)-selectedRange(1);327   aLut = aLut/(selectedRange(2)-selectedRange(1));328 else329   % remap from 0..1 to 0..numBins-1330   if numBins ~= 1331     binStep = 1/(numBins-1);332     start = ceil(selectedRange(1)/binStep);333     stop  = floor(selectedRange(2)/binStep);334     k = start+1:stop+1;335     aLut(k) = 0:1/(length(k)-1):1;336   else337     aLut(1) = 0; %in case someone specifies numBins = 1, which is just silly338   end339 end340 341 imgTileRow=1;342 for k=1:numTiles(1)+1343   if k == 1  %special case: top row344     imgTileNumRows = dimTile(1)/2; %always divisible by 2 because of padding345     mapTileRows = [1 1];346   else 347     if k == numTiles(1)+1 %special case: bottom row      348       imgTileNumRows = dimTile(1)/2;349       mapTileRows = [numTiles(1) numTiles(1)];350     else %default values351       imgTileNumRows = dimTile(1); 352       mapTileRows = [k-1, k]; %[upperRow lowerRow]353     end354   end355   356   % loop over columns of the tileMappings cell array357   imgTileCol=1;358   for l=1:numTiles(2)+1359     if l == 1 %special case: left column360       imgTileNumCols = dimTile(2)/2;361       mapTileCols = [1, 1];362     else363       if l == numTiles(2)+1 % special case: right column364         imgTileNumCols = dimTile(2)/2;365         mapTileCols = [numTiles(2), numTiles(2)];366       else %default values367         imgTileNumCols = dimTile(2);368         mapTileCols = [l-1, l]; % right left369       end370     end371     372     % Extract four tile mappings373     ulMapTile = tileMappings{mapTileRows(1), mapTileCols(1)};374     urMapTile = tileMappings{mapTileRows(1), mapTileCols(2)};375     blMapTile = tileMappings{mapTileRows(2), mapTileCols(1)};376     brMapTile = tileMappings{mapTileRows(2), mapTileCols(2)};377 378     % Calculate the new greylevel assignments of pixels 379     % within a submatrix of the image specified by imgTileIdx. This 380     % is done by a bilinear interpolation between four different mappings 381     % in order to eliminate boundary artifacts.382     383     normFactor = imgTileNumRows*imgTileNumCols; %normalization factor  384     imgTileIdx = {imgTileRow:imgTileRow+imgTileNumRows-1, ...385                  imgTileCol:imgTileCol+imgTileNumCols-1};386 387     imgPixVals = grayxform(I(imgTileIdx{1},imgTileIdx{2}), aLut);388     389     % calculate the weights used for linear interpolation between the390     % four mappings391     rowW = repmat((0:imgTileNumRows-1)‘,1,imgTileNumCols);392     colW = repmat(0:imgTileNumCols-1,imgTileNumRows,1);393     rowRevW = repmat((imgTileNumRows:-1:1)‘,1,imgTileNumCols);394     colRevW = repmat(imgTileNumCols:-1:1,imgTileNumRows,1);395     396     claheI(imgTileIdx{1}, imgTileIdx{2}) = ...397         (rowRevW .* (colRevW .* double(grayxform(imgPixVals,ulMapTile)) + ...398                      colW    .* double(grayxform(imgPixVals,urMapTile)))+ ...399          rowW    .* (colRevW .* double(grayxform(imgPixVals,blMapTile)) + ...400                      colW    .* double(grayxform(imgPixVals,brMapTile))))...401         /normFactor;402     403     imgTileCol = imgTileCol + imgTileNumCols;    404   end %over tile cols405   imgTileRow = imgTileRow + imgTileNumRows;406 end %over tile rows407 408 %-----------------------------------------------------------------------------409 410 function [I, selectedRange, fullRange, numTiles, dimTile, clipLimit,...411           numBins, noPadRect, distribution, alpha, ...412           int16ClassChange] = parseInputs(varargin)413 414 narginchk(1,13);415 416 I = varargin{1};417 validateattributes(I, {‘uint8‘, ‘uint16‘, ‘double‘, ‘int16‘, ‘single‘}, ...418               {‘real‘, ‘2d‘, ‘nonsparse‘, ‘nonempty‘}, ...419               mfilename, ‘I‘, 1);420 421 % convert int16 to uint16422 if isa(I,‘int16‘)423   I = int16touint16(I);424   int16ClassChange = true;425 else426   int16ClassChange = false;427 end428 429 if any(size(I) < 2)430   error(message(‘images:adapthisteq:inputImageTooSmall‘))431 end432 433 %Other options434 %%%%%%%%%%%%%%435 436 %Set the defaults437 distribution = ‘uniform‘;438 alpha   = 0.4;439 440 if isa(I, ‘double‘) || isa(I,‘single‘)441   fullRange = [0 1];442 else443   fullRange(1) = I(1);         %copy class of the input image444   fullRange(1:2) = [-Inf Inf]; %will be clipped to min and max445   fullRange = double(fullRange);446 end447 448 selectedRange   = fullRange;449 450 %Set the default to 256 bins regardless of the data type;451 %the user can override this value at any time452 numBins = 256;453 normClipLimit = 0.01;454 numTiles = [8 8];455 456 checkAlpha = false;457 458 validStrings = {‘NumTiles‘,‘ClipLimit‘,‘NBins‘,‘Distribution‘,...459                 ‘Alpha‘,‘Range‘};460 461 if nargin > 1462   done = false;463 464   idx = 2;465   while ~done466     input = varargin{idx};467     inputStr = validatestring(input, validStrings,mfilename,‘PARAM‘,idx);468 469     idx = idx+1; %advance index to point to the VAL portion of the input 470 471     if idx > nargin472       error(message(‘images:adapthisteq:missingValue‘, inputStr))473     end474     475     switch inputStr476 477      case ‘NumTiles‘478        numTiles = varargin{idx};479        validateattributes(numTiles, {‘double‘}, {‘real‘, ‘vector‘, ...480                            ‘integer‘, ‘finite‘,‘positive‘},...481                      mfilename, inputStr, idx);482 483        if (any(size(numTiles) ~= [1,2]))484          error(message(‘images:adapthisteq:invalidNumTilesVector‘, inputStr))485        end486        487        if any(numTiles < 2)488          error(message(‘images:adapthisteq:invalidNumTilesValue‘, inputStr))489        end490       491      case ‘ClipLimit‘492       normClipLimit = varargin{idx};493       validateattributes(normClipLimit, {‘double‘}, ...494                     {‘scalar‘,‘real‘,‘nonnegative‘},...495                     mfilename, inputStr, idx);496       497       if normClipLimit > 1498         error(message(‘images:adapthisteq:invalidClipLimit‘, inputStr))499       end500      501      case ‘NBins‘502       numBins = varargin{idx};      503       validateattributes(numBins, {‘double‘}, {‘scalar‘,‘real‘,‘integer‘,...504                           ‘positive‘}, mfilename, ‘NBins‘, idx);505      506      case ‘Distribution‘507       validDist = {‘rayleigh‘,‘exponential‘,‘uniform‘};508       distribution = validatestring(varargin{idx}, validDist, mfilename,...509                                   ‘Distribution‘, idx);510      511      case ‘Alpha‘512       alpha = varargin{idx};513       validateattributes(alpha, {‘double‘},{‘scalar‘,‘real‘,...514                           ‘nonnan‘,‘positive‘,‘finite‘},...515                     mfilename, ‘Alpha‘,idx);516       checkAlpha = true;517 518      case ‘Range‘519       validRangeStrings = {‘original‘,‘full‘};520       rangeStr = validatestring(varargin{idx}, validRangeStrings,mfilename,...521                               ‘Range‘,idx);522       523       if strmatch(rangeStr,‘original‘)524         selectedRange = double([min(I(:)), max(I(:))]);525       end526      527      otherwise528       error(message(‘images:adapthisteq:inputString‘)) %should never get here529     end530     531     if idx >= nargin532       done = true;533     end534     535     idx=idx+1;536   end537 end538 539 540 %% Pre-process the inputs541 %%%%%%%%%%%%%%%%%%%%%%%%%%542 543 dimI = size(I);544 dimTile = dimI ./ numTiles;545 546 %check if tile size is reasonable547 if any(dimTile < 1)548   error(message(‘images:adapthisteq:inputImageTooSmallToSplit‘, num2str( numTiles )))549 end550 551 if checkAlpha552   if strcmp(distribution,‘uniform‘)553     error(message(‘images:adapthisteq:alphaShouldNotBeSpecified‘, distribution))554   end555 end556 557 %check if the image needs to be padded; pad if necessary;558 %padding occurs if any dimension of a single tile is an odd number559 %and/or when image dimensions are not divisible by the selected 560 %number of tiles561 rowDiv  = mod(dimI(1),numTiles(1)) == 0;562 colDiv  = mod(dimI(2),numTiles(2)) == 0;563 564 if rowDiv && colDiv565   rowEven = mod(dimTile(1),2) == 0;566   colEven = mod(dimTile(2),2) == 0;  567 end568 569 noPadRect = [];570 if  ~(rowDiv && colDiv && rowEven && colEven)571   padRow = 0;572   padCol = 0;573   574   if ~rowDiv575     rowTileDim = floor(dimI(1)/numTiles(1)) + 1;576     padRow = rowTileDim*numTiles(1) - dimI(1);577   else578     rowTileDim = dimI(1)/numTiles(1);579   end580   581   if ~colDiv582     colTileDim = floor(dimI(2)/numTiles(2)) + 1;583     padCol = colTileDim*numTiles(2) - dimI(2);584   else585     colTileDim = dimI(2)/numTiles(2);586   end587   588   %check if tile dimensions are even numbers589   rowEven = mod(rowTileDim,2) == 0;590   colEven = mod(colTileDim,2) == 0;591   592   if ~rowEven593     padRow = padRow+numTiles(1);594   end595 596   if ~colEven597     padCol = padCol+numTiles(2);598   end599   600   padRowPre  = floor(padRow/2);601   padRowPost = ceil(padRow/2);602   padColPre  = floor(padCol/2);603   padColPost = ceil(padCol/2);604   605   I = padarray(I,[padRowPre  padColPre ],‘symmetric‘,‘pre‘);606   I = padarray(I,[padRowPost padColPost],‘symmetric‘,‘post‘);607 608   %UL corner (Row, Col), LR corner (Row, Col)609   noPadRect.ulRow = padRowPre+1;610   noPadRect.ulCol = padColPre+1;611   noPadRect.lrRow = padRowPre+dimI(1);612   noPadRect.lrCol = padColPre+dimI(2);613 end614 615 %redefine this variable to include the padding616 dimI = size(I);617 618 %size of the single tile619 dimTile = dimI ./ numTiles;620 621 %compute actual clip limit from the normalized value entered by the user622 %maximum value of normClipLimit=1 results in standard AHE, i.e. no clipping;623 %the minimum value minClipLimit would uniformly distribute the image pixels624 %across the entire histogram, which would result in the lowest possible625 %contrast value626 numPixInTile = prod(dimTile);627 minClipLimit = ceil(numPixInTile/numBins);628 clipLimit = minClipLimit + round(normClipLimit*(numPixInTile-minClipLimit));629 630 %-----------------------------------------------------------------------------