Histogram Normalization implementation

[Creator]

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;
}


[Creator]

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:

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

Equalized Histogram:(the histogram is almost constant)

[Creator]

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

[Creator]

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();
}


[Creator]

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;
   }


[Creator]

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;
   }


[Administrator]

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