Entropy thresholding

//Uses concept of Disorder of image (Entropy)
//Calculates disorder of image for each treshold level, then picks minimum disorder
//Author: m.sugrue @ rhul.ac.uk, 2004

ImageSpace ImageSpace::EntropyThreshold(void)
{ 
ImageSpace *ReturnImageSpace = new ImageSpace;
BYTE *b1 = new BYTE[ XMax*YMax ];
 
int ThresholdLevel, *hist; 
double prob[256], Sum_prob_1k = 0, Sum_prob_kl = 0, 
  Sum_prob_ln_1k = 0, Sum_prob_ln_kl = 0, Entropy[256];

 hist = Histogram();    //gets histogram of image 

 for( int i = 0; i< 256; i++)
  prob[i] =  (double) hist[i] / (XMax * YMax);  // probability distribution
 

 for( int k = 0; k< 255; k++)  // k is trial threshold level
 {
  Sum_prob_1k = 0; Sum_prob_kl = 0; 
  Sum_prob_ln_1k = 0; Sum_prob_ln_kl = 0;
  
  for(i = 1; i < k; i++)   //From 1 to k
  {        //Sums probablities to k
   Sum_prob_1k += prob[i];  //Sums prob times Log of prob to k
   if(prob[i]!=0) Sum_prob_ln_1k += (prob[i] * log(prob[i]));
  }
  
  for(i = k; i < 256; i++) //From k to end
  {        //Sums prob of k to end
   Sum_prob_kl += prob[i];  //Sums prob times log of prob
   if(prob[i]!=0) Sum_prob_ln_kl += (prob[i] * log(prob[i]));
  }
     //Final equation of entropy for each k
  Entropy[k] = log(Sum_prob_1k) + log(Sum_prob_kl)
   - (Sum_prob_ln_1k / Sum_prob_1k)
   - (Sum_prob_ln_kl / Sum_prob_kl);
     //protects against divide by zero
  if(Entropy[k]<0) Entropy[k]=0;
 }  

 ThresholdLevel = 0;
 for(k = 0; k<256; k++) //Finds Maximum
  if(Entropy[k] > Entropy[ThresholdLevel]) ThresholdLevel = k;
       //Thresholds there
  cout << "Thresholded at :" << ThresholdLevel << endl;
 for( int y=0; y < YMax; y++ ) for( int x=0; x < XMax; x++ ) 
  b1[ x + y * XMax ] = ( b0[ x + y * XMax ] > ThresholdLevel ) ? 255 : 0;
 
 ReturnImageSpace->Load(b1);
 IsThreshold = TRUE;

return *ReturnImageSpace;
}

//Here is the same code in c#
//A 'Frame' is simply a grid of pixel values (float) and label values (int)

  public override Bitmap ShowMotion()
  {
   //Entropy threshold

   Frame frame = new Frame(_MotionFrame[_Decomposition - 1].ToBitmap());

   int ThresholdLevel;
   double[] prob = new double[256], Entropy = new double[256];
   double Sum_prob_1k = 0, Sum_prob_kl = 0,
    Sum_prob_ln_1k = 0, Sum_prob_ln_kl = 0;

   //get Histogram
   int[] hist = new int[256];
   for(int i = 0; i < frame.Width; i++)
    for(int j = 0; j < frame.Height; j++)
     hist[ (int) frame[i, j].Intensity ]++;

   //probability distribution
   double normalisationfactor = (frame.Width * frame.Height);
   for(int i = 0; i < 256; i++)
    prob[i] = (double) hist[ i ] / normalisationfactor;

   for( int k = 0; k< 255; k++)  // k is trial threshold level
   {
    Sum_prob_1k = 0; Sum_prob_kl = 0; 
    Sum_prob_ln_1k = 0; Sum_prob_ln_kl = 0;
  
    for(int i = 1; i < k; i++)   //From 1 to k
    {        //Sums probablities to k
     Sum_prob_1k += prob[i];  //Sums prob times Log of prob to k
     if(prob[i]!=0) Sum_prob_ln_1k += (prob[i] * Math.Log(prob[i]));
    }
  
    for(int i = k; i < 256; i++)  //From k to end
    {        //Sums prob of k to end
     Sum_prob_kl += prob[i];  //Sums prob times log of prob
     if(prob[i]!=0) Sum_prob_ln_kl += (prob[i] * Math.Log(prob[i]));
    }
    //Final equation of entropy for each k
    Entropy[k] = Math.Log(Sum_prob_1k) + Math.Log(Sum_prob_kl)
     - (Sum_prob_ln_1k / Sum_prob_1k)
     - (Sum_prob_ln_kl / Sum_prob_kl);
    //protects against divide by zero
    if(Entropy[k]<0) Entropy[k]=0;
   }  

   ThresholdLevel = 0;
   for(int k = 0; k<256; k++) //Finds Maximum
    if(Entropy[k] > Entropy[ThresholdLevel]) ThresholdLevel = k;

   //Thresholds there
   Trace.WriteLine("Thresholded at :" + ThresholdLevel.ToString() );

   for(int i=0; i < frame.Width; i++ ) 
    for(int j=0; j < frame.Height; j++ ) 
     frame[i, j].Label = ( frame[i, j].Intensity  > ThresholdLevel) ? 0 : -1;

   return frame.Labels_ToBitmap();
  }