Make an Adaptive Histogram Equalization

[snibgo]

My purpose for filling nulls is to compensate for the fact that my camera doesn't generate 16-bit sRGB values. It generates either 8-bit sRGB (in jpg files) or 14-bit RGB raw files, which I then convert with dcraw and IM to create 16-bit tiff files. The histograms of these files will have some counts of zero, associated with an artificially high count at one or both ends of the stream of zeros. The RGB data might have nulls spread evenly throughout the pixel intensities, but the RGB->sRGB conversion will redistribute nulls.

As an example, part of the histogram might have counts "25, 25, 25, 25" where another has values "0, 0, 0, 100". A naive search for peaks will find a peak at "100" although it probably isn't a real peak.

To counteract this problem, I identify a stream of (n) (where n > 0) entries that have zero counts, followed by a count of at least (n+1). I redistribute the entire count over the (n+1) entries. So "0, 0, 0, 100" would be redistributed as "25, 25, 25, 25".

(If the histogram were an image, some form of blurring would be equivalent.)

In the code, ph->Values[x].Count is the histogram count for value x where 0 <= x < 65536.

Code: Select all

  int nFills = 0;
  for (i = 0; i < NUM_VALUES-1; i++) {
    ValueT * pvi = &ph->Values[i];
    if (!pvi->Count) {
      // Find next non-zero and test if its count is more than distance.
      for (j = i+1; j < NUM_VALUES; j++) {
        ValueT * pvj = &ph->Values[j];
        if (pvj->Count) {
          int n = j - i + 1;
          if (pvj->Count >= n) {
            nFills++;
            int count = pvj->Count / n;
            int rem = pvj->Count - (n * count);
            //printf ("Back filling %i * %i from %i to %i, remainder %i\n", n, count, j, i, rem);
            // We also spread remainder.
            int kLim = j - rem + 1;
            for (k = i; k <= j; k++) {
              ValueT * pvk = &ph->Values[k];
              if (k < j && pvk->Count) printf ("FillNullCounts: ** %g\n", pvk->Count);
              pvk->Count = count;
              if (k >= kLim) pvk->Count++;
              //printf ("= %i\n", pvk->Count);
            }

            i = j;
          }
          break; // out of j loop
        }
      }
    }
  }

There are some issues with this.

1. The code is grossly slow when I have large blocks of zero counts.

2. I should probably redistribute even where count < n+1.

3. My code redistributes remainders, so "0, 0, 0, 106" will become "26, 26, 27, 27". However, my count array (ph->Values[x].Count) is now floating-point (to cater for alpha masking), so this should change to give "26.5, 26.5, 26.5, 26.5".

4. I always redistribute downwards. So a stream of "0, 0, 0, 100, 0, 0, 0, 4" will become "25, 25, 25, 25, 1, 1, 1, 1", which is probably not ideal.

5. A consequence of issue 4 is that a count at histogram[0] is never redistributed, so there is often an incorrect peak there. If the first entries are "78, 0, 0, 0, 0, 100", they will become "78, 25, 25, 25".

[fmw42]

That is quite an elaborate method. You are essentially spreading out a value to fill a bunch of preceding zeros. Thus keeping the total count the same. Interesting.

[snibgo]

It needs to be fairly elaborate, so I can identify peaks and troughs in the histogram. Here's an example where the image has been automatically partitioned into four areas, corresponding to the four peaks in the histogram of the blurred image. Each segment spans the tones from one histogram trough to the next.



Identifying peaks and troughs is quite complex. My current method uses a moving window across the histogram counts. The window width is 65536 * 0.1, extending each side of the histogram entry under question. A count is a peak if three conditions are true:

- this count is at least as large than any other in the window
- this count is at least mean_count_in_window * 1.1
- this count is at least mean_count_in_window + 2

I then ignore small peaks, defined as having a count less than the highest count * 0.1, although this wasn't necessary in this image.

A trough is the lowest value between two peaks.

[fmw42]

That is very clever. You might be interesting in http://stackoverflow.com/questions/3790 ... -detection. I used it to code a peak finder for my unperspective script for finding the corners of the picture in polar coordinates.

[snibgo]

Thanks for the link. My first stab at the problem came up with something very similar; we know a peak has occurred when we have dropped by a certain delta, and a trough has occured after we have risen by a certain delta. It's very fast, needing just a single pass through the data.

The trouble was, I couldn't find a value for delta that worked for different images. It would find two peaks very close together that was really just one peak with noise, or would miss a small but clear peak in a flattish histogram, or would falsely claim peaks in a noisy histogram.

I realised that, for my purpose, what distinguished two peaks wasn't the depth of the trough between them, but the horizontal distance between them. Doubtless a hybrid approach is possible, perhaps with automatic setting of depth or window-width.

[fmw42]

It is a tough problem. I had to add tests for spacing in width and height to get rid of some false peaks and also throw out small peaks.

Glad you found an approach that worked for you.

[snibgo]

The technique works for me but, in fairness, I can't claim it is perfect. If an entry qualifies as a peak, I record it and resume searching half-a-window further on. Many other entries within the window often have the same count, but I ignore them. It might be better to choose the most central of all the equal counts; the difficulty is that they each have their own different windows, so might not qualify as a peak.

Furthermore, if we imagine a weird histogram that has a constant low count then suddenly rises to a constant high count and remains like that to the end, my algorithm will claim that the first high count is a peak. Fortunately, real-world photos don't have such weird histograms.

[fmw42]

I have now folded in my redist script a contrast gain limiting technique similar to that used by user snibgo in his equalize code. Thanks to user snibgo for all the useful comments and explanations.