Histogram matching of two Images without using histeq

Question

It is well known that histeq in MATLAB can perform histogram matching so that an image's histogram is transformed to look like another histogram. I am trying to perform this same operation without using histeq. I'm aware that you need to calculate the CDFs between the two images, but I'm not sure what to do next. What do I do?

Answer 1

Histogram matching is concerned with transforming one image's histogram so that it looks like another. The basic principle is to compute the histogram of each image individually, then compute their discrete cumulative distribution functions (CDFs). Let's denote the CDF of first image as while the CDF of the second image is . Therefore, would denote what the CDF value is for intensity x for the first image.

Once you calculate the CDFs for each of the images, you need to compute a mapping that transforms one intensity from the first image so that it is in agreement with the intensity distribution of the second image. To do this, for each intensity in the first image - let's call this which will be from [0,255] assuming an 8-bit image - we must find an intensity in the second image (also in the range of [0,255]) such that:

There may be a case where we won't get exactly an equality, so what you would need to do is find the smallest absolute difference between and . In other words, for a mapping M, for each entry of , we must find an intensity such that:

You would do this for all 256 values, and we would produce a mapping. Once you find this mapping, you simply have to apply this mapping on the first image to get it to look like the intensity distribution of the second image. A rough (and perhaps inefficient) algorithm would look something like this. Let im1 be the first image (of type uint8) while im2 is the second image (of type uint8):

M = zeros(256,1,'uint8'); %// Store mapping - Cast to uint8 to respect data type
hist1 = imhist(im1); %// Compute histograms
hist2 = imhist(im2);
cdf1 = cumsum(hist1) / numel(im1); %// Compute CDFs
cdf2 = cumsum(hist2) / numel(im2);

%// Compute the mapping
for idx = 1 : 256
    [~,ind] = min(abs(cdf1(idx) - cdf2));
    M(idx) = ind-1;
end

%// Now apply the mapping to get first image to make
%// the image look like the distribution of the second image
out = M(double(im1)+1);

out should contain your matched image where it transforms the intensity distribution of the first image to match that of the second image. Take special care of the out statement. The intensity range of im1 spans between [0,255], but MATLAB's indexing for arrays starts at 1. Therefore, we need to add 1 to every value of im1 so we can properly index into M to produce our output. However, im1 is of type uint8, and MATLAB saturates values should you try and go beyond 255. As such, to ensure that we get to 256, we must cast to a data type that is beyond 8-bit precision. I decided to use double, then when we add 1 to every value in im1, we will span between 1 to 256 so we can properly index into M. Also take not that when I find the location that minimizes the difference, I also must subtract by 1 as the data type spans from [0,255].