Contrast Limited Adaptive Histogram Equalization with Gamma distribution

Hello all,

I would like to implement an extension of the well-known CLAHE histogram equalization method in Matlab. Currently, Matlab provides routines to locally equalize the histogram of the image while mapping it to:

 1. A uniform distribution
 2. A Rayleigh distribution
 3. An exponential distribution

By digging into the code of adapthisteq (the function for CLAHE in matlab) I've found the following pieces of code that specify which distribution to map the output histogram:

      case 'uniform',
      scale =  valSpread/numPixInTile;
      mapping = min(selectedRange(1) + histSum*scale,...
                  selectedRange(2)); %limit to max

      case 'rayleigh', % suitable for underwater imagery
      % pdf = (t./alpha^2).*exp(-t.^2/(2*alpha^2))*U(t)
      % cdf = 1-exp(-t.^2./(2*alpha^2))
      hconst = 2*alpha^2;
      vmax = 1 - exp(-1/hconst);
      val = vmax*(histSum/numPixInTile);
      val(val>=1) = 1-eps; % avoid log(0)
      temp = sqrt(-hconst*log(1-val));
      mapping = min(selectedRange(1)+temp*valSpread,...
                  selectedRange(2)); %limit to max

       case 'exponential',
       % pdf = alpha*exp(-alpha*t)*U(t)
       % cdf = 1-exp(-alpha*t)
       vmax = 1 - exp(-alpha);
       val = (vmax*histSum/numPixInTile);
       val(val>=1) = 1-eps;
       temp = -1/alpha*log(1-val);
       mapping = min(selectedRange(1)+temp*valSpread,selectedRange(2));

This makes me think that temp is just the inverse function of the corresponding cdf when t=1, i.e., of vmax. Then, I guess I should just find the inverse of the gamma cummulative distribution, and set it in temp, and that should do the job. However, I am not being able to work out such an inverse function. The gamma probability and cummulative distributions are defined as:

pdf(t; k, theta) = (t^{k-1} exp(-t/theta)) / (theta^k gamma(k))

cdf(t; k, theta) = gammainc(k, 1/theta) / gamma(k)

where gamma is the gamma function, and gammainc is the incomplete gamma function, both defined on matlab. Any clue about how to solve this?

Thank you very much!

Answer 1

[Image Analyst]

Correct, to equalize you do use inverse transform sampling where the CDF is computed. Some functions have analytical forms for the inverse cdf, and if so, you can use those, like I did for my attached function that draws from the Rayleigh distribution.

%-----------------------------------------------------------------
% KEY PART, RIGHT HERE!!!!
% Invert the CDF of the Rayleigh function to get a function that can
% generate random numbers drawn from a Rayleigh distribution,
% given numbers drawn from a uniform distribution.
rayleighDistNumbers = sqrt(-log(1-uniformlyDistributedRandomNumbers)*(2*sigma^2));
%-----------------------------------------------------------------

They're doing the same thing as I do - there is a known formula for the inverse. I just took the formula as a given. To invert it yourself analytically would probably involve some heavy duty math that I didn't want to bother doing when I had the answer already in hand.

[Adrián]

Thanks for your quick response. I have found that there is a built-in function in matlab that allows one to compute the inverse of the cumulative gamma distribution:

http://es.mathworks.com/help/stats/gaminv.html

I am trying with that, and see if I can get the equalization I am looking for.

Thanks!