Selective histogram equalization (only on a specified area of the image)

Question

I am developping on Qt creator, with opencv.

I have to developp a program that does the histogram equalization of an image. My images are 16bits grayscale images so I cannot use the opencv function "equalizeHist" because it only works with 8bit grayscale images.

The code I wrote to do it is the following :

void equalizeHist_16U(Mat* img)
{
    long hist[65535] = {0};
    double ratio;
    int i, j;

    assert(img->channels() == 1);
    assert(img->type() == CV_16U);

    ratio = 65535.0 / (img->cols*img->rows);

    //Cumulative histogram calculation
    compute_hist_16U(img, hist, true);

    for(i=0 ; i<img->cols ; i++)
    {
        for(j=0 ; j<img->rows ; j++)
        {
            img->at<unsigned short>(j,i) = ratio*hist[img->at<unsigned short>(j,i)];
        }
    }

}

long compute_hist_16U (Mat* img, long* hist, bool cumul)
{
    unsigned short i, j, k;
    long* b;
    long max = 0;

    //is the image 16bits grayscale ?
    assert(img->channels() == 1);
    assert(CV_16U == img->type());

    //histogram calculation
    for(i=0 ; i<img->cols ; i++)
    {
        for(j=0 ; j<img->rows ; j++)
        {
            hist[img->at<unsigned short>(j,i)]++;
            if(hist[img->at<unsigned short>(j,i)] > max)
                max = hist[img->at<unsigned short>(j,i)];
        }
    }

    //Cumulative histogram calculation (if cumul=true)
    if(cumul)
    {
        for(b=hist ; b<hist+65535 ; b++)
        {
            *(b+1) += *b;
        }
    }
    return (cumul ? hist[65535] : max);
}

It does what I expected, now I want to do the histogram equialization of my image but only on a specified part of the image. I added x1,x2,y1,y2 parameters to my function and changed the bounds of my "for" like this (the lines of code i changed have arrows):

---->void equalizeHist_16U(Mat* img, int x1, int x2, int y1, int y2)
{
    long hist[65535] = {0};
    double ratio;
    int i, j;

    assert(img->channels() == 1);
    assert(img->type() == CV_16U);

    ratio = 65535.0 / (img->cols*img->rows);

    //Cumulative histogram calculation
    compute_hist_16U(img, hist, true);

    ---->for(i=x1 ; i<=x2 ; i++)
    {
        ---->for(j=y1 ; j<=y2 ; j++)
        {
            img->at<unsigned short>(j,i) = ratio*hist[img->at<unsigned short>(j,i)];
        }
    }

}

---->long compute_hist_16U (Mat* img, long* hist, bool cumul, int x1, int x2, int y1, int y2)
{
    unsigned short i, j, k;
    long* b;
    long max = 0;

    //is the image 16bits grayscale ?
    assert(img->channels() == 1);
    assert(CV_16U == img->type());

    //histogram calculation
    ---->for(i=x1 ; i<=x2  ; i++)
    {
        ---->for(j=y1 ; j<=y2 ; j++)
        {
            hist[img->at<unsigned short>(j,i)]++;
            if(hist[img->at<unsigned short>(j,i)] > max)
                max = hist[img->at<unsigned short>(j,i)];
        }
    }

    //Cumulative histogram calculation (if cumul=true)
    if(cumul)
    {
        for(b=hist ; b<hist+65535 ; b++)
        {
            *(b+1) += *b;
        }
    }
    return (cumul ? hist[65535] : max);
}

But it does not work as expected, my image is not equalized, I do not have extrem values (clear white and dark black) on my image. If I try

equalizeHist_16U(&img, 0, 50, 0, 50)

The image i get is very very bright And if I try

equalizeHist(&img, 300, 319, 220, 239)

The image I get is very very dark

I think I have made a mistake in the loops bounds but I can't find where ! Maybe you have an idea ?

Thank you in advance

Antonio Antonio · Accepted Answer · 2015-05-07T22:12:56

Preliminary:
Did you notice you are not using at all your second version of the cumulative histogram function?

void equalizeHist_16U(Mat* img, int x1, int x2, int y1, int y2)

is calling

compute_hist_16U(img, hist, true);

and not:

long compute_hist_16U (Mat* img, long* hist, bool cumul, int x1, int x2, int y1, int y2)

(I suppose you wanted to post the last one, otherwise I wouldn't see why you posted the code :) )

Actual answer:

Everything can become much easier if you use cv::Mat rois, through the operator ().

Your function would become as following:

void equalizeHist_16U(Mat* img, int x1, int x2, int y1, int y2) {

   //Here you should check you have x2 > x1 and y2 > y1 and y1,x1 >0 and x2 <= img->width and y2 <= img->height

   cv::Rect roi(x1,y1,x2-x1,y2-y1); 
   //To reproduce exactly the behaviour you seem to target,
   //it should be x2-x2+1 and y2-y1+1. 
   //But you should get used on the fact that extremes are,
   //as a convention, excluded

   cv::Mat& temp = *img; //Otherwise using operator() is complicated
   cv::Mat roiMat = temp(roi); //This doesn't do any memory copy, just creates a new header!
   void equalizeHist_16U(&roiMat); //Your original function!!

}

And that's it! If this doesn't work, then it means your original function that process the whole image has a bug that you couldn't notice before.

When I'll have a bit of time I'll post a couple of suggestions to make your function faster (e.g., you should avoid using .at, you should compute the max value in your histogram at the end of the histogram computation, you should create a look up table of short where you multiply your histogram by ratio, so that applying the histogram becomes much faster; and instead of the ratio variable which causes floating point conversion, you could simply divide your histogram elements by the constant (img->width*img->height)) and more neat (you should pass Mat by reference, not using a pointer, that's C-style, not C++)

Furthermore:

Why do you return a value from compute_hist_16U?
long hist[65535] should long hist[65536], so that index 65535 is valid. First of all, 65535 is the white value in your image. Furthermore, you use it in your cycle, when you have b+1 with b=hist+65534 (the last cycle)
```
for(b=hist ; b<hist+65535 ; b++)
{
    *(b+1) += *b;
}
```

Selective histogram equalization (only on a specified area of the image)

2 Answers