local histogram equalization in matlab / python

Question

i am very new in matlab. i want to write the code for local histogram equalization . i have been written code for global histogram equalization and i know that local equalization means do equalization for each part of image seperately but my question is that how i should choose this part of images ? for example should i do equalization for each 100 pixel that are neighbor separate of other pixels ? in the other word how i can take apart image to some parts and then do equalization to each part?

Answer 1

The most naive way to do what you ask is split up your image into non-overlapping blocks, do your global histogram code on that block and save it to the output. Suppose you defined the rows and columns of these non-overlapping blocks as the variables rows and cols. In your case, let's say it's 100 x 100, so rows = 100; cols = 100;. You would simply loop over each non-overlapping block, do your histogram equalization then set this to the same locations in the output.

Something like this below, assuming your image is stored in im:

rows = 100;
cols = 100;
out = zeros(size(im)); % Declare output variable

for ii = 1 : rows : size(im, 1)
    for jj = 1 : cols : size(im, 2)
        % Get the block
        row_begin = ii;
        row_end = min(size(im, 1), ii + rows);
        col_begin = jj;
        col_end = min(size(im, 2), jj + cols);
        blk = im(row_begin : row_end, col_begin : col_end, :);
        % Perform histogram equalization with the block stored in blk
        % ...

        % Assume the output of this is stored in O
        out(row_begin : row_end, col_begin : col_end, :) = O;
    end
end

Note the intricacy of the variable blk that stores the non-overlapping block. We let the beginning row and column simply be the loop counter ii and jj, but the ending row and column we must make sure that it's bounded by the dimensions of the image. That's why the min call is there. Otherwise, the ending row and column is simply the beginning row and column added by the size of the block in the corresponding dimensions. Also note that I've used : to index into the third dimension in case you have a colour image. Grayscale should not affect this code. You finally need to use the same indexing when storing the output in the output image. Note that I've assumed this is stored in the variable O which is the output of your customized histogram equalization function.

The output out will contain your locally histogram equalized image. Take note that you could theoretically do this in one line using blockproc in the image processing toolbox if you have it. This processes distinct blocks in your image and applies some function to it. Assuming your histogram equalization function is called hsteq, you would simply do this:

rows = 100; cols = 100;
out = blockproc(im, [rows, cols], @(s) hsteq(s.data));

The first input is the image you want to process, the second input defines the block size and finally the last element is the function you want to apply to each block. Note that blockproc supplies a customized structure into your function and so what is important is that you pull out the data field in the structure. This should produce the same output as the code above with loops.

Answer 2

We can use the tile-based local (adaptive) histogram equalization to implement AHE (as suggested in the other answer), but in that case we need to implement a bilinear interpolation-like technique to prevent sudden change of contrasts at the edges of the window, e.g., observe the equalized output below with python implementation of the same (here a 50x50 window is used for the tile):

def AHE(im, tile_x=8, tile_y=8):        
    h, w = im.shape
    out = np.zeros(im.shape) # Declare output variable

    for i in range(0, h, tile_x):
        for j in range(0, w, tile_y):
            # Get the block
            blk = im[i: min(i + tile_x, h), j: min(j + tile_y, w)]
            probs = get_distr(blk)
            out[i: min(i + tile_x, h), j: min(j + tile_y, w)] = CHE(blk, probs)
    return out

def CHE(im, probs):
    T = np.array(list(map(int, 255*np.cumsum(probs))))
    return T[im]

def get_distr(im):
   hist, _ = np.histogram(im.flatten(),256,[0,256])
   return hist / hist.sum()

We could instead implement the AHE algorithm from this thesis:

The implementation of algorithm yields better results (without the boundary artifacts):