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

/*

 * ANSI C code from the article

 * "Contrast Limited Adaptive Histogram Equalization"

 * by Karel Zuiderveld, karel@cv.ruu.nl

 * in "Graphics Gems IV", Academic Press, 1994

 *

 *

 *  These functions implement Contrast Limited Adaptive Histogram Equalization.

 *  The main routine (CLAHE) expects an input image that is stored contiguously in

 *  memory;  the CLAHE output image overwrites the original input image and has the

 *  same minimum and maximum values (which must be provided by the user).

 *  This implementation assumes that the X- and Y image resolutions are an integer

 *  multiple of the X- and Y sizes of the contextual regions. A check on various other

 *  error conditions is performed.

 *

 *  #define the symbol BYTE_IMAGE to make this implementation suitable for

 *  8-bit images. The maximum number of contextual regions can be redefined

 *  by changing uiMAX_REG_X and/or uiMAX_REG_Y; the use of more than 256

 *  contextual regions is not recommended.

 *

 *  The code is ANSI-C and is also C++ compliant.

 *

 *  Author: Karel Zuiderveld, Computer Vision Research Group,

 *         Utrecht, The Netherlands (karel@cv.ruu.nl)

 */

#ifdef BYTE_IMAGE

typedef unsigned char kz_pixel_t;     /* for 8 bit-per-pixel images */

#define uiNR_OF_GREY (256)

#else

typedef unsigned short kz_pixel_t;     /* for 12 bit-per-pixel images (default) */

# define uiNR_OF_GREY ()

#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 = ;      /* max. # contextual regions in x-direction */

const unsigned int uiMAX_REG_Y = ;      /* 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 -;       /* # of regions x-direction too large */

    if (uiNrY > uiMAX_REG_Y) return -;       /* # of regions y-direction too large */

    if (uiXRes % uiNrX) return -;      /* x-resolution no multiple of uiNrX */

    if (uiYRes & uiNrY) return -;      /* y-resolution no multiple of uiNrY */

    if (Max >= uiNR_OF_GREY) return -;       /* maximum too large */

    if (Min >= Max) return -;          /* minimum equal or larger than maximum */

    if (uiNrX <  || uiNrY < ) return -;/* at least 4 contextual regions required */

    if (fCliplimit == 1.0) return ;      /* is OK, immediately returns original image. */

    if (uiNrBins == ) uiNrBins = ;      /* default value when not specified */

    pulMapArray=(unsigned long *)malloc(sizeof(unsigned long)*uiNrX*uiNrY*uiNrBins);

    if (pulMapArray == ) return -;      /* 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<<;          /* 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 = , pImPointer = pImage; uiY < uiNrY; uiY++) {

    for (uiX = ; 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 - ) * uiXRes;          /* skip lines, set pointer */

    }

    /* Interpolate greylevel mappings to get CLAHE image */

    for (pImPointer = pImage, uiY = ; uiY <= uiNrY; uiY++) {

    if (uiY == ) {                      /* special case: top row */

        uiSubY = uiYSize >> ;  uiYU = ; uiYB = ;

    }

    else {

        if (uiY == uiNrY) {                  /* special case: bottom row */

        uiSubY = uiYSize >> ;    uiYU = uiNrY-;     uiYB = uiYU;

        }

        else {                      /* default values */

        uiSubY = uiYSize; uiYU = uiY - ; uiYB = uiYU + ;

        }

    }

    for (uiX = ; uiX <= uiNrX; uiX++) {

        if (uiX == ) {                  /* special case: left column */

        uiSubX = uiXSize >> ; uiXL = ; uiXR = ;

        }

        else {

        if (uiX == uiNrX) {              /* special case: right column */

            uiSubX = uiXSize >> ;  uiXL = uiNrX - ; uiXR = uiXL;

        }

        else {                      /* default values */

            uiSubX = uiXSize; uiXL = uiX - ; uiXR = uiXL + ;

        }

        }

        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 - ) * uiXRes;

    }

    free(pulMapArray);                      /* free space for histograms */

    return ;                          /* 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 = ;  pulBinPointer = pulHistogram;

    for (i = ; i < uiNrGreylevels; i++) { /* calculate total number of excess pixels */

    lBinExcess = (long) pulBinPointer[i] - (long) ulClipLimit;

    if (lBinExcess > ) 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 = ; 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 < ) ulStepSize = ;          /* 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 = ; i < uiNrGreylevels; i++) pulHistogram[i] = 0L; /* clear histogram */

    for (i = ; 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 = ;

    const float fScale = ((float)(Max - Min)) / ulNrOfPixels;

    const unsigned long ulMin = (unsigned long) Min;

    for (i = ; 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) ( + (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 = ;

    if (uiNum & (uiNum - ))   /* If uiNum is not a power of two, use division */

    for (uiYCoef = , uiYInvCoef = uiYSize; uiYCoef < uiYSize;

     uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {

    for (uiXCoef = , 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 >>= ) uiShift++;           /* Calculate 2log of uiNum */

    for (uiYCoef = , uiYInvCoef = uiYSize; uiYCoef < uiYSize;

         uiYCoef++, uiYInvCoef--,pImage+=uiIncr) {

         for (uiXCoef = , 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);

        }

    }

    }

}

 if (pulHistogram[i] > ulUpper)

             {       /* high bin count */

                 ulNrExcess -= (ulClipLimit -pulHistogram[i]); pulHistogram[i]=ulClipLimit;

             }   

 function out = adapthisteq(varargin)

 %ADAPTHISTEQ Contrast-limited Adaptive Histogram Equalization (CLAHE).

 %   ADAPTHISTEQ enhances the contrast of images by transforming the

 %   values in the intensity image I.  Unlike HISTEQ, it operates on small

 %   data regions (tiles), rather than the entire image. Each tile's

 %   contrast is enhanced, so that the histogram of the output region

 %   approximately matches the specified histogram. The neighboring tiles

 %   are then combined using bilinear interpolation in order to eliminate

 %   artificially induced boundaries.  The contrast, especially

 %   in homogeneous areas, can be limited in order to avoid amplifying the

 %   noise which might be present in the image.

 %

 %   J = ADAPTHISTEQ(I) Performs CLAHE on the intensity image I.

 %

 %   J = ADAPTHISTEQ(I,PARAM1,VAL1,PARAM2,VAL2...) sets various parameters.

 %   Parameter names can be abbreviated, and case does not matter. Each

 %   string parameter is followed by a value as indicated below:

 %

 %   'NumTiles'     Two-element vector of positive integers: [M N].

 %                  [M N] specifies the number of tile rows and

 %                  columns.  Both M and N must be at least .

 %                  The total number of image tiles is equal to M*N.

 %

 %                  Default: [ ].

 %

 %   'ClipLimit'    Real scalar from  to .

 %                  'ClipLimit' limits contrast enhancement. Higher numbers

 %                  result in more contrast.

 %

 %                  Default: 0.01.

 %

 %   'NBins'        Positive integer scalar.

 %                  Sets number of bins for the histogram used in building a

 %                  contrast enhancing transformation. Higher values result

 %                  in greater dynamic range at the cost of slower processing

 %                  speed.

 %

 %                  Default: .

 %

 %   'Range'        One of the strings: 'original' or 'full'.

 %                  Controls the range of the output image data. If 'Range'

 %                  is set to 'original', the range is limited to

 %                  [min(I(:)) max(I(:))]. Otherwise, by default, or when

 %                  'Range' is set to 'full', the full range of the output

 %                  image class is used (e.g. [ ] for uint8).

 %

 %                  Default: 'full'.

 %

 %   'Distribution' Distribution can be one of three strings: 'uniform',

 %                  'rayleigh', 'exponential'.

 %                  Sets desired histogram shape for the image tiles, by

 %                  specifying a distribution type.

 %

 %                  Default: 'uniform'.

 %

 %   'Alpha'        Nonnegative real scalar.

 %                  'Alpha' is a distribution parameter, which can be supplied

 %                  when 'Dist' is set to either 'rayleigh' or 'exponential'.

 %

 %                  Default: 0.4.

 %

 %   Notes

 %   -----

 %   - 'NumTiles' specify the number of rectangular contextual regions (tiles)

 %     into which the image is divided. The contrast transform function is

 %     calculated for each of these regions individually. The optimal number of

 %     tiles depends on the type of the input image, and it is best determined

 %     through experimentation.

 %

 %   - The 'ClipLimit' is a contrast factor that prevents over-saturation of the

 %     image specifically in homogeneous areas.  These areas are characterized

 %     by a high peak in the histogram of the particular image tile due to many

 %     pixels falling inside the same gray level range. Without the clip limit,

 %     the adaptive histogram equalization technique could produce results that,

 %     in some cases, are worse than the original image.

 %

 %   - ADAPTHISTEQ can use Uniform, Rayleigh, or Exponential distribution as

 %     the basis for creating the contrast transform function. The distribution

 %     that should be used depends on the type of the input image.

 %     For example, underwater imagery appears to look more natural when the

 %     Rayleigh distribution is used.

 %

 %   Class Support

 %   -------------

 %   Intensity image I can be uint8, uint16, int16, double, or single.

 %   The output image J has the same class as I.

 %

 %   Example

 %   ---------

 %   Apply Contrast-Limited Adaptive Histogram Equalization to an

 %   image and display the results.

 %

 %      I = imread('tire.tif');

 %      A = adapthisteq(I,'clipLimit',0.02,'Distribution','rayleigh');

 %      figure, imshow(I);

 %      figure, imshow(A);

 %

 %   Example

 %   ---------

 %

 %   Apply Contrast-Limited Adaptive Histogram Equalization to a color

 %   photograph.

 %

 %      [X MAP] = imread('shadow.tif');

 %      RGB = ind2rgb(X,MAP); % convert indexed image to truecolor format

 %      cform2lab = makecform('srgb2lab');

 %      LAB = applycform(RGB, cform2lab); %convert image to L*a*b color space

 %      L = LAB(:,:,)/; % scale the values to range from  to

 %      LAB(:,:,) = adapthisteq(L,'NumTiles',[ ],'ClipLimit',0.005)*;

 %      cform2srgb = makecform('lab2srgb');

 %      J = applycform(LAB, cform2srgb); %convert back to RGB

 %      figure, imshow(RGB); %display the results

 %      figure, imshow(J);

 %

 %   See also HISTEQ.

 %   Copyright - The MathWorks, Inc.

 %   $Revision: 1.1.6.12 $  $Date: // :: $

 %   References:

 %      Karel Zuiderveld, "Contrast Limited Adaptive Histogram Equalization",

 %      Graphics Gems IV, p. -, code: p. -

 %

 %      Hanumant Singh, Woods Hole Oceanographic Institution, personal

 %      communication

 %--------------------------- The algorithm ----------------------------------

 %

 %  . Obtain all the inputs:

 %    * image

 %    * number of regions in row and column directions

 %    * number of bins for the histograms used in building image transform

 %      function (dynamic range)

 %    * clip limit for contrast limiting (normalized from  to )

 %    * other miscellaneous options

 %  . Pre-process the inputs:

 %    * determine real clip limit from the normalized value

 %    * if necessary, pad the image before splitting it into regions

 %  . Process each contextual region (tile) thus producing gray level mappings

 %    * extract a single image region

 %    * make a histogram for this region using the specified number of bins

 %    * clip the histogram using clip limit

 %    * create a mapping (transformation function) for this region

 %  . Interpolate gray level mappings in order to assemble final CLAHE image

 %    * extract cluster of four neighboring mapping functions

 %    * process image region partly overlapping each of the mapping tiles

 %    * extract a single pixel, apply four mappings to that pixel, and

 %      interpolate between the results to obtain the output pixel; repeat

 %      over the entire image

 %

 %  See code for further details.

 %

 %-----------------------------------------------------------------------------

 [I, selectedRange, fullRange, numTiles, dimTile, clipLimit, numBins, ...

  noPadRect, distribution, alpha, int16ClassChange] = parseInputs(varargin{:});

 tileMappings = makeTileMappings(I, numTiles, dimTile, numBins, clipLimit, ...

                                 selectedRange, fullRange, distribution, alpha);

 %Synthesize the output image based on the individual tile mappings.

 out = makeClaheImage(I, tileMappings, numTiles, selectedRange, numBins,...

                      dimTile);

 if int16ClassChange

   % Change uint16 back to int16 so output has same class as input.

   out = uint16toint16(out);

 end

 if ~isempty(noPadRect) %do we need to remove padding?

   out = out(noPadRect.ulRow:noPadRect.lrRow, ...

             noPadRect.ulCol:noPadRect.lrCol);

 end

 %-----------------------------------------------------------------------------

 function tileMappings = ...

     makeTileMappings(I, numTiles, dimTile, numBins, clipLimit,...

                      selectedRange, fullRange, distribution, alpha)

 numPixInTile = prod(dimTile);

 tileMappings = cell(numTiles);

 % extract and process each tile

 imgCol = ;

 for col=:numTiles(),

   imgRow = ;

   for row=:numTiles(),

     tile = I(imgRow:imgRow+dimTile()-,imgCol:imgCol+dimTile()-);

     % for speed, call MEX file directly thus avoiding costly

     % input parsing of imhist

     tileHist = imhistc(tile, numBins, , fullRange()); 

     tileHist = clipHistogram(tileHist, clipLimit, numBins);

     tileMapping = makeMapping(tileHist, selectedRange, fullRange, ...

                               numPixInTile, distribution, alpha);

     % assemble individual tile mappings by storing them in a cell array;

     tileMappings{row,col} = tileMapping;

     imgRow = imgRow + dimTile();

   end

   imgCol = imgCol + dimTile(); % move to the next column of tiles

 end

 %-----------------------------------------------------------------------------

 % Calculate the equalized lookup table (mapping) based on cumulating the input

 % histogram.  Note: lookup table is rescaled in the selectedRange [Min..Max].

 function mapping = makeMapping(imgHist, selectedRange, fullRange, ...

                                numPixInTile, distribution, alpha)

 histSum = cumsum(imgHist);

 valSpread  = selectedRange() - selectedRange();

 switch distribution

  case 'uniform',

   scale =  valSpread/numPixInTile;

   mapping = min(selectedRange() + histSum*scale,...

                 selectedRange()); %limit to max

  case 'rayleigh', % suitable for underwater imagery

   % pdf = (t./alpha^).*exp(-t.^/(*alpha^))*U(t)

   % cdf = -exp(-t.^./(*alpha^))

   hconst = *alpha^;

   vmax =  - exp(-/hconst);

   val = vmax*(histSum/numPixInTile);

   val(val>=) = -eps; % avoid log()

   temp = sqrt(-hconst*log(-val));

   mapping = min(selectedRange()+temp*valSpread,...

                 selectedRange()); %limit to max

  case 'exponential',

   % pdf = alpha*exp(-alpha*t)*U(t)

   % cdf = -exp(-alpha*t)

   vmax =  - exp(-alpha);

   val = (vmax*histSum/numPixInTile);

   val(val>=) = -eps;

   temp = -/alpha*log(-val);

   mapping = min(selectedRange()+temp*valSpread,selectedRange());

  otherwise,

   error(message('images:adapthisteq:distributionType')) %should never get here

 end

 %rescale the result to be between  and  for later use by the GRAYXFORM

 %private mex function

 mapping = mapping/fullRange();

 %-----------------------------------------------------------------------------

 % This function clips the histogram according to the clipLimit and

 % redistributes clipped pixels across bins below the clipLimit

 function imgHist = clipHistogram(imgHist, clipLimit, numBins)

 % total number of pixels overflowing clip limit in each bin

 totalExcess = sum(max(imgHist - clipLimit,));  

 % clip the histogram and redistribute the excess pixels in each bin

 avgBinIncr = floor(totalExcess/numBins);

 upperLimit = clipLimit - avgBinIncr; % bins larger than this will be

                                      % set to clipLimit

 % this loop should speed up the operation by putting multiple pixels

 % into the "obvious" places first

 for k=:numBins

   if imgHist(k) > clipLimit

     imgHist(k) = clipLimit;

   else

     if imgHist(k) > upperLimit % high bin count

       totalExcess = totalExcess - (clipLimit - imgHist(k));

       imgHist(k) = clipLimit;

     else

       totalExcess = totalExcess - avgBinIncr;

       imgHist(k) = imgHist(k) + avgBinIncr;

     end

   end

 end

 % this loops redistributes the remaining pixels, one pixel at a time

 k = ;

 while (totalExcess ~= )

   %keep increasing the step as fewer and fewer pixels remain for

   %the redistribution (spread them evenly)

   stepSize = max(floor(numBins/totalExcess),);

   for m=k:stepSize:numBins

     if imgHist(m) < clipLimit

       imgHist(m) = imgHist(m)+;

       totalExcess = totalExcess - ; %reduce excess

       if totalExcess ==

         break;

       end

     end

   end

   k = k+; %prevent from always placing the pixels in bin #

   if k > numBins % start over if numBins was reached

     k = ;

   end

 end

 %-----------------------------------------------------------------------------

 % This function interpolates between neighboring tile mappings to produce a

 % new mapping in order to remove artificially induced tile borders.

 % Otherwise, these borders would become quite visible.  The resulting

 % mapping is applied to the input image thus producing a CLAHE processed

 % image.

 function claheI = makeClaheImage(I, tileMappings, numTiles, selectedRange,...

                                  numBins, dimTile)

 %initialize the output image to zeros (preserve the class of the input image)

 claheI = I;

 claheI(:) = ;

 %compute the LUT for looking up original image values in the tile mappings,

 %which we created earlier

 if ~(isa(I,'double') || isa(I,'single'))

   k = selectedRange()+ : selectedRange()+;

   aLut = zeros(length(k),);

   aLut(k) = (k-)-selectedRange();

   aLut = aLut/(selectedRange()-selectedRange());

 else

   % remap from .. to ..numBins-

   if numBins ~=

     binStep = /(numBins-);

     start = ceil(selectedRange()/binStep);

     stop  = floor(selectedRange()/binStep);

     k = start+:stop+;

     aLut(k) = :/(length(k)-):;

   else

     aLut() = ; %in case someone specifies numBins = , which is just silly

   end

 end

 imgTileRow=;

 for k=:numTiles()+

   if k ==   %special case: top row

     imgTileNumRows = dimTile()/; %always divisible by  because of padding

     mapTileRows = [ ];

   else

     if k == numTiles()+ %special case: bottom row

       imgTileNumRows = dimTile()/;

       mapTileRows = [numTiles() numTiles()];

     else %default values

       imgTileNumRows = dimTile();

       mapTileRows = [k-, k]; %[upperRow lowerRow]

     end

   end

   % loop over columns of the tileMappings cell array

   imgTileCol=;

   for l=:numTiles()+

     if l ==  %special case: left column

       imgTileNumCols = dimTile()/;

       mapTileCols = [, ];

     else

       if l == numTiles()+ % special case: right column

         imgTileNumCols = dimTile()/;

         mapTileCols = [numTiles(), numTiles()];

       else %default values

         imgTileNumCols = dimTile();

         mapTileCols = [l-, l]; % right left

       end

     end

     % Extract four tile mappings

     ulMapTile = tileMappings{mapTileRows(), mapTileCols()};

     urMapTile = tileMappings{mapTileRows(), mapTileCols()};

     blMapTile = tileMappings{mapTileRows(), mapTileCols()};

     brMapTile = tileMappings{mapTileRows(), mapTileCols()};

     % Calculate the new greylevel assignments of pixels

     % within a submatrix of the image specified by imgTileIdx. This

     % is done by a bilinear interpolation between four different mappings

     % in order to eliminate boundary artifacts.

     normFactor = imgTileNumRows*imgTileNumCols; %normalization factor

     imgTileIdx = {imgTileRow:imgTileRow+imgTileNumRows-, ...

                  imgTileCol:imgTileCol+imgTileNumCols-};

     imgPixVals = grayxform(I(imgTileIdx{},imgTileIdx{}), aLut);

     % calculate the weights used for linear interpolation between the

     % four mappings

     rowW = repmat((:imgTileNumRows-)',1,imgTileNumCols);

     colW = repmat(:imgTileNumCols-,imgTileNumRows,);

     rowRevW = repmat((imgTileNumRows:-:)',1,imgTileNumCols);

     colRevW = repmat(imgTileNumCols:-:,imgTileNumRows,);

     claheI(imgTileIdx{}, imgTileIdx{}) = ...

         (rowRevW .* (colRevW .* double(grayxform(imgPixVals,ulMapTile)) + ...

                      colW    .* double(grayxform(imgPixVals,urMapTile)))+ ...

          rowW    .* (colRevW .* double(grayxform(imgPixVals,blMapTile)) + ...

                      colW    .* double(grayxform(imgPixVals,brMapTile))))...

         /normFactor;

     imgTileCol = imgTileCol + imgTileNumCols;

   end %over tile cols

   imgTileRow = imgTileRow + imgTileNumRows;

 end %over tile rows

 %-----------------------------------------------------------------------------

 function [I, selectedRange, fullRange, numTiles, dimTile, clipLimit,...

           numBins, noPadRect, distribution, alpha, ...

           int16ClassChange] = parseInputs(varargin)

 narginchk(,);

 I = varargin{};

 validateattributes(I, {'uint8', 'uint16', 'double', 'int16', 'single'}, ...

               {'real', '2d', 'nonsparse', 'nonempty'}, ...

               mfilename, 'I', );

 % convert int16 to uint16

 if isa(I,'int16')

   I = int16touint16(I);

   int16ClassChange = true;

 else

   int16ClassChange = false;

 end

 if any(size(I) < )

   error(message('images:adapthisteq:inputImageTooSmall'))

 end

 %Other options

 %%%%%%%%%%%%%%

 %Set the defaults

 distribution = 'uniform';

 alpha   = 0.4;

 if isa(I, 'double') || isa(I,'single')

   fullRange = [ ];

 else

   fullRange() = I();         %copy class of the input image

   fullRange(:) = [-Inf Inf]; %will be clipped to min and max

   fullRange = double(fullRange);

 end

 selectedRange   = fullRange;

 %Set the default to  bins regardless of the data type;

 %the user can override this value at any time

 numBins = ;

 normClipLimit = 0.01;

 numTiles = [ ];

 checkAlpha = false;

 validStrings = {'NumTiles','ClipLimit','NBins','Distribution',...

                 'Alpha','Range'};

 if nargin >

   done = false;

   idx = ;

   while ~done

     input = varargin{idx};

     inputStr = validatestring(input, validStrings,mfilename,'PARAM',idx);

     idx = idx+; %advance index to point to the VAL portion of the input 

     if idx > nargin

       error(message('images:adapthisteq:missingValue', inputStr))

     end

     switch inputStr

      case 'NumTiles'

        numTiles = varargin{idx};

        validateattributes(numTiles, {'double'}, {'real', 'vector', ...

                            'integer', 'finite','positive'},...

                      mfilename, inputStr, idx);

        if (any(size(numTiles) ~= [,]))

          error(message('images:adapthisteq:invalidNumTilesVector', inputStr))

        end

        if any(numTiles < )

          error(message('images:adapthisteq:invalidNumTilesValue', inputStr))

        end

      case 'ClipLimit'

       normClipLimit = varargin{idx};

       validateattributes(normClipLimit, {'double'}, ...

                     {'scalar','real','nonnegative'},...

                     mfilename, inputStr, idx);

       if normClipLimit >

         error(message('images:adapthisteq:invalidClipLimit', inputStr))

       end

      case 'NBins'

       numBins = varargin{idx};

       validateattributes(numBins, {'double'}, {'scalar','real','integer',...

                           'positive'}, mfilename, 'NBins', idx);

      case 'Distribution'

       validDist = {'rayleigh','exponential','uniform'};

       distribution = validatestring(varargin{idx}, validDist, mfilename,...

                                   'Distribution', idx);

      case 'Alpha'

       alpha = varargin{idx};

       validateattributes(alpha, {'double'},{'scalar','real',...

                           'nonnan','positive','finite'},...

                     mfilename, 'Alpha',idx);

       checkAlpha = true;

      case 'Range'

       validRangeStrings = {'original','full'};

       rangeStr = validatestring(varargin{idx}, validRangeStrings,mfilename,...

                               'Range',idx);

       if strmatch(rangeStr,'original')

         selectedRange = double([min(I(:)), max(I(:))]);

       end

      otherwise

       error(message('images:adapthisteq:inputString')) %should never get here

     end

     if idx >= nargin

       done = true;

     end

     idx=idx+;

   end

 end

 %% Pre-process the inputs

 %%%%%%%%%%%%%%%%%%%%%%%%%%

 dimI = size(I);

 dimTile = dimI ./ numTiles;

 %check if tile size is reasonable

 if any(dimTile < )

   error(message('images:adapthisteq:inputImageTooSmallToSplit', num2str( numTiles )))

 end

 if checkAlpha

   if strcmp(distribution,'uniform')

     error(message('images:adapthisteq:alphaShouldNotBeSpecified', distribution))

   end

 end

 %check if the image needs to be padded; pad if necessary;

 %padding occurs if any dimension of a single tile is an odd number

 %and/or when image dimensions are not divisible by the selected

 %number of tiles

 rowDiv  = mod(dimI(),numTiles()) == ;

 colDiv  = mod(dimI(),numTiles()) == ;

 if rowDiv && colDiv

   rowEven = mod(dimTile(),) == ;

   colEven = mod(dimTile(),) == ;

 end

 noPadRect = [];

 if  ~(rowDiv && colDiv && rowEven && colEven)

   padRow = ;

   padCol = ;

   if ~rowDiv

     rowTileDim = floor(dimI()/numTiles()) + ;

     padRow = rowTileDim*numTiles() - dimI();

   else

     rowTileDim = dimI()/numTiles();

   end

   if ~colDiv

     colTileDim = floor(dimI()/numTiles()) + ;

     padCol = colTileDim*numTiles() - dimI();

   else

     colTileDim = dimI()/numTiles();

   end

   %check if tile dimensions are even numbers

   rowEven = mod(rowTileDim,) == ;

   colEven = mod(colTileDim,) == ;

   if ~rowEven

     padRow = padRow+numTiles();

   end

   if ~colEven

     padCol = padCol+numTiles();

   end

   padRowPre  = floor(padRow/);

   padRowPost = ceil(padRow/);

   padColPre  = floor(padCol/);

   padColPost = ceil(padCol/);

   I = padarray(I,[padRowPre  padColPre ],'symmetric','pre');

   I = padarray(I,[padRowPost padColPost],'symmetric','post');

   %UL corner (Row, Col), LR corner (Row, Col)

   noPadRect.ulRow = padRowPre+;

   noPadRect.ulCol = padColPre+;

   noPadRect.lrRow = padRowPre+dimI();

   noPadRect.lrCol = padColPre+dimI();

 end

 %redefine this variable to include the padding

 dimI = size(I);

 %size of the single tile

 dimTile = dimI ./ numTiles;

 %compute actual clip limit from the normalized value entered by the user

 %maximum value of normClipLimit= results in standard AHE, i.e. no clipping;

 %the minimum value minClipLimit would uniformly distribute the image pixels

 %across the entire histogram, which would result in the lowest possible

 %contrast value

 numPixInTile = prod(dimTile);

 minClipLimit = ceil(numPixInTile/numBins);

 clipLimit = minClipLimit + round(normClipLimit*(numPixInTile-minClipLimit));

 %-----------------------------------------------------------------------------