Histogram Normalization implementation

Microsoft Visual C++
User avatar
Creator
Posts: 157
Joined: Tue Dec 16, 2008, 20:52
Location: Hannover, Germany
Contact:

Histogram Normalization implementation

Postby Creator » Sat Apr 23, 2011, 18:18

We have Histogram Equalization in Graff px implemented.
It is appeared that in Internet almost impossible to find a Histogram Normalization algorithm. So, here I propose a C++ implementation which uses naive OpenCV data containers. The algorithm is straight-forward but should be easily understandable. I suggest to make a small competition on optimizing it :-)

Code: Select all

void Normalization()
{
   register int   x, y;
   register int   i, j;
   unsigned int   nC[256], n = 0;
   double         pC[256], pN[256];
   double         cC[256], cN[256];

   for (i=0; i < 256; i++) {
      nC[i] = 0;
      cC[i] = 0;
   }
   n = this->cols * this->rows;

   //calculate the nubber of occurence
   vector<Mat> channels(this->channels());
   split(*this, channels);
   for (x = 0; x < this->cols; x++)
      for (y = 0; y < this->rows; y++) {
         int val = 0;
         for (i = 0; i < this->channels(); i++)   val += channels[i].at<unsigned char>(y, x);
         val /= this->channels();
         nC[(unsigned char) val]++;
      }

   //calculating probabilities
   for (i=0; i < 256; i++) {
      pC[i] = (double) nC[i]/ (double) n;
      pN[i] = exp(-pow((double) i - 127, 2) / (2*32*32)) / sqrt(2*Pi*32*32);
      //calculatinc cumulative distribution
      for (j=0; j<=i; j++) {
         cC[i] += pC[j];
         cN[i] += pN[j];
      }
   }

   //for (i=0; i < 256; i++) cN[i] = (double ) i / 256;   //uncomment this to achieve histogram equalization
   
   for (x = 0; x < this->cols; x++)
      for (y = 0; y < this->rows; y++)
         for (i = 0; i < this->channels(); i++) {
            double arg = cC[channels[i].at<unsigned char>(y,x)];
            double Min = abs(cN[0] - arg);
            int MinJ = 0;
            for (j = 1; j < 256; j++)
               if (Min > abs(cN[j] - arg)) {
                  Min = abs(cN[j] - arg);
                  MinJ = j;
               }
            channels[i].at<unsigned char>(y, x) = (unsigned char) MinJ;
         }
               
   merge(channels, *this);
   channels.clear();
   return 0;
}

User avatar
Creator
Posts: 157
Joined: Tue Dec 16, 2008, 20:52
Location: Hannover, Germany
Contact:

Re: Histogram Normalization implementation

Postby Creator » Wed May 04, 2011, 20:51

Some explanations to what the specification of the algorithm is. The function modifies the histogram of some descreete signal (a digital image in out case). After the modidication the histogram should obey the normal gaussian distribution law.
To show that on example, please concider the images with its histogrms in the right lower corners
Original image:
cat_org.jpg

Stretched Histogram: (please note, that the form of the histogram keeps the same)
cat_scl.jpg

Equalized Histogram:(the histogram is almost constant)
cat_qul.jpg

User avatar
Creator
Posts: 157
Joined: Tue Dec 16, 2008, 20:52
Location: Hannover, Germany
Contact:

Re: Histogram Normalization implementation

Postby Creator » Wed May 04, 2011, 20:52

Normalized Histogram:(the histogram looks like gaussian "bell" function)
cat_nor.jpg

User avatar
Creator
Posts: 157
Joined: Tue Dec 16, 2008, 20:52
Location: Hannover, Germany
Contact:

Re: Histogram Normalization implementation

Postby Creator » Wed May 04, 2011, 22:14

I suggest to split the code in few parts to make the task easier. Moreover, we can add the look-up tables for optimisation:
1. Variables Declaration and Initialization

Code: Select all

void Normalization()
{
   register int x, y;
   register int i, j;
   int mu, sigma;                  // Gaussian parameters   unsigned int lut[256];            // Look-Up Table (LUT)
   unsigned int nC[256], n;         // Histogram of "this" image, number of pixels in image
   double pC[256], cC[256];         // Image intencity value occurance probability, Cumalative Distribution Function of pC
   double pN[256], cN[256];         // Normal intencity value occurance probability, Cumalative Distribution Function of pN
   
   memset(nC, 0, 256*sizeof(unsigned int));
   n = this->cols * this->rows;
   mu = 256 / 2; sigma = mu / 3;


2. Calculating the Histogram:

Code: Select all

   std::vector<Mat> vChannels(this->channels());
   split(*this, vChannels);
   for (x = 0; x < this->cols; x++)
      for (y = 0; y < this->rows; y++) {
         int val = 0;
         for (i = 0; i < this->channels(); i++) val += vChannels[i].at<unsigned char>(y, x);
         val /= this->channels();   // Graff px style
         nC[val]++;
      }


3. Calculating the Probabilities:

Code: Select all

   for (i=0; i < 256; i++) {
      pC[i] = static_cast<double>(nC[i]) / n;
      // Calculating cumulative distribution
      cC[i] = (i==0) ? pC[i] : cC[i-1] + pC[i];
   }

   // Creating the form of the normal histogram, and corresponding CDF
   for(i = 0; i < 256; i++) {
      pN[i] = exp(-pow((double) i - mu, 2) / (2*sigma*sigma)) / sqrt(2*Pi*sigma*sigma);
      // Calculating cumulative distribution
      cN[i] = (i==0) ? pN[i] : cN[i-1] + pN[i];
   }


4.Creating look-up table (LUT):

Code: Select all

   for (i = 0; i < 256; i++) {
      double Min = abs(cN[0] - cC[i]);
      int MinJ = 0;
      for (j = 1; j < 256; j++)
         if (Min > abs(cN[j] - cC[i])) {
            Min = abs(cN[j] - cC[i]);
            MinJ = j;
         }
      lut[i] = MinJ;
   }
   


5. Resulting Image Transforming:

Code: Select all

 
   for (x = 0; x < this->cols; x++)
      for (y = 0; y < this->rows; y++)
         for (i = 0; i < this->channels(); i++) {
            unsigned char val = vChannels[i].at<unsigned char>(y, x);
            vChannels[i].at<unsigned char>(y, x) = static_cast<unsigned char>(lut[val]);
         }
   merge(vChannels, *this);       
   vChannels.clear();
}

User avatar
Creator
Posts: 157
Joined: Tue Dec 16, 2008, 20:52
Location: Hannover, Germany
Contact:

Re: Histogram Normalization implementation

Postby Creator » Tue May 10, 2011, 16:12

To optimize №4, it was suggested to use the binary search algorithm. Its complexity is O(log n). Its implementation is here:

4.Creating look-up table (LUT):

Code: Select all

   for (i = 0; i < 256; i++) {
      int first = 0;      
      int last =255; 
      int mid;
      while ( first < last ) {
         mid = first + (last-first) / 2;
         if (cC[i] <= cN[mid]) last = mid;
         else  first = mid + 1;
      }
       lut[i] = last;
   }

User avatar
Creator
Posts: 157
Joined: Tue Dec 16, 2008, 20:52
Location: Hannover, Germany
Contact:

Re: Histogram Normalization implementation

Postby Creator » Tue May 10, 2011, 16:19

Also, thanks to Rotsor the algorithm was furter optimized to the complexoty O(n):

4.Creating look-up table (LUT):

Code: Select all

   j = 0;
   for (i = 0; i < 256; i++) {
      while((j < 255) && (abs(cN[j + 1] - cC[i]) <= abs(cN[j] - cC[i]))) j++;
      lut[i] = j;
   }
Last edited by Creator on Wed May 18, 2011, 19:50, edited 1 time in total.

Administrator
Site Admin
Posts: 20
Joined: Tue Dec 16, 2008, 20:29
Contact:

Re: Histogram Normalization implementation

Postby Administrator » Fri May 13, 2011, 14:52

For those who prefere evaluation the algorithm in milliseconds instead of notations of its complexity, I've evaluated the algorithm for look-up table estimation (the code part №4). Evealuation was performed in a Intel Core i7 950 Processor 3,07GHz with 24GB DDR3 RAM. The test image had the resolution 571x700 pixels with 3 color bands, coded with 16 bits each (so the length of the histogram is not 256 as coded above, but 65 536).
  • Agorithm O(n^2): 32,8 ms
  • Agorithm O(n log n): 2,2 ms
  • Agorithm O(n): 0,06 ms
So, the finally the speed-up appeared to be more than 500 times :!:


Return to “Programming”

Who is online

Users browsing this forum: No registered users and 13 guests