//From http://www.richelbilderbeek.nl/CppDoHistogramEqualization.htm
void DoHistogramEqualization(const TImage * const source, TImage * const target)
{
  assert(source!=0 && "Source image is NULL");
  assert(target!=0 && "Target image is NULL");
  assert(source->Picture->Bitmap!=0 && "Source bitmap is NULL");
  assert(target->Picture->Bitmap!=0 && "Target bitmap is NULL");
  assert(source->Picture->Bitmap->PixelFormat == pf24bit && "Source bitmap must be 24 bit");
  assert(target->Picture->Bitmap->PixelFormat == pf24bit && "Target bitmap must be 24 bit");
  //Get the width and height from the source
  const int width  = source->Picture->Bitmap->Width;
  const int height = source->Picture->Bitmap->Height;
  //Set the target's width and height
  target->Picture->Bitmap->Width  = width;
  target->Picture->Bitmap->Height = height;

  const int surface = width * height;
  const int nGreyValues = 256; //There are 256 different pixel intensities
  const std::vector<int> histogram = GetImageHistogram(source);
  assert(nGreyValues==static_cast<int>(histogram.size()));
  const std::vector<int> cumulativeHistogram = GetCumulativeHistogram(histogram);
  assert(nGreyValues==static_cast<int>(cumulativeHistogram.size()));

  //Works, but anybody knows how to use std::for_each or std::transform for this?
  std::vector<int> rescaledHistogram(nGreyValues,0);
  for (int i=0; i!=nGreyValues; ++i)
  {
    //'surface + 1' to prevent that rescaledGreyValue == 256
    const int rescaledGreyValue
      = static_cast<int>(
          static_cast<double>(nGreyValues)
        * static_cast<double>(cumulativeHistogram[i])
        / static_cast<double>(surface + 1) );
    assert(rescaledGreyValue >= 0);
    assert(rescaledGreyValue < 256);
    rescaledHistogram[i] = rescaledGreyValue;
  }

  for (int y=0; y!=height; ++y)
  {
    const unsigned char * lineSource
      = static_cast<const unsigned char *>(
        source->Picture->Bitmap->ScanLine[y]);
    unsigned char * lineTarget
      = static_cast<unsigned char *>(
        target->Picture->Bitmap->ScanLine[y]);
    for (int x=0; x!=width; ++x)
    {
      const int greyOriginal
        = (lineSource[x*3+0] + lineSource[x*3+1] + lineSource[x*3+2]) / 3;
      assert(greyOriginal >=   0);
      assert(greyOriginal  < 256);
      const int greyNew = rescaledHistogram[greyOriginal];
      assert(greyNew >= 0);
      assert(greyNew  < 256);

      lineTarget[x*3+0]=greyNew; //Blue
      lineTarget[x*3+1]=greyNew; //Green
      lineTarget[x*3+2]=greyNew; //Red
    }
  }
}