# Histogram matching (specification) in Python

I'm trying to implement an algorithm in which I first pad each row of the image with a fixed amount of new pixels in a certain range, apply Gaussian smoothing to the row cumulative histograms in vertical (y) direction, and thus obtain new cumulative histograms for each row in the end.

After having obtained the new row cumulative histograms, my task is to get back/restore an intensity image from the cumulative histograms. Specifically, I would like to match each row histogram to the new, corresponding row histogram (to those which I obtained after padding and Gaussian filtering). Therefore, if there are, for example, 1000 rows in an image I will match 1000 row histograms and restore the original image.

My code in Python is below, where I came until the step where I need to perform histogram matching but I got stuck there:

import cv2
import numpy as np
import scipy.ndimage as ndi
from matplotlib import pyplot as plt
from scipy.interpolate import interp1d

# Threshold of intensity values for padding
T = 200
rows,cols = img.shape
# Number of pixels to pad each bin of a row histogram
N = 30
########################################################################
for i in range(0,rows):
img_row = img[i,]
# Calculate the row histogram (without padding)
# Calculate the cumulative row histogram (without padding)
# Copy the row to prepare for padding
# Apply uniform padding to all the bins less than T
# Calculate the cumulative row histogram of the padded row
# Accumulate the cumulative histograms of padded rows in a matrix

# Apply 1D-Gaussian filtering on the padded cumulative histogram along the columns
Gauss_cdf = ndi.gaussian_filter1d(cdf_hist_Padded, sigma=2, axis=0, output=np.float64, mode='nearest')

# Normalize all the CDFs to get values between [0,1]
norm_Gauss_cdf = Gauss_cdf/Gauss_cdf.max()

# Take the first original and padded+smoothed row cumulative histogram
H_hat_zero = norm_Gauss_cdf[0,:]

# I would like to match now 'H_zero' and 'H_hat_zero'. When I find how to do that,
# I  will apply this to all rows in a loop. How do I perform the matching?..


If I understood right, you are stuck in matching a given histogram into a desired one and creating a new image from this matched histogram obtained by your filtering method.

I would first suggest you to get rid of all the unnecessary stuff (including the python code) above and isolate your problem as "histogram matching" in mathematical terms.

Now I will describe how you can approximately match a given histogram to a desired one in two setps: by first converting it into that of a uniform (equalized) one, and in then converting this uniform one to the desired one. Finally we shall combine these two steps to get the answer.

The method is based on converting a random variable into another by means of a transform $G(\cdot)$. You will replace the random variable $X$ with image intensity function $I(n,m)$.

Consider a uniform random variable $X$ whose CDF is $F_X(x) = U(x) = P(X < x) = x$. We want to generate a new rv $Y$ from $X$, by means of a transform $Y=G(X)$, such that its CDF is the desired one. We shall therefore find this transform $G(x)$.

Under suitable conditions: $$F_Y(y)= P(Y< y) = P(G(X)<y) = P(X< G^{-1}(y)) = F_X(G^{-1}(y))=G^{-1}(y)$$

from which we deduce that the transform that converts uniform $X$ to arbitrary $Y$ with a desired CDF is $G(x) = F_Y^{-1}(x)$ from which we obtain: $$\boxed{Y = F_Y^{-1}(X)}$$

Therefore the transform that converts a uniform RV $X$ to an arbitrary (desired) RV $Y$ is $$\boxed{G(x) = F_Y^{-1}(x)}$$

The transform that converts an arbitrary RV $Y$ to a uniform RV $X$ is $$\boxed{ G^{-1}(x)=F_Y(x) }$$

Finaly, the transform that converts arbitrary RV $Y$ to another arbitrary RV $Z=T(Y)$ is $$\boxed{ T(y)=F_Z^{-1}(F_Y(y)) }$$

We are almost done, we will finish by converting above algorithm to Image manipulation and also clarify how you will overcome the problem of not having those CDF's as formulas.

Replace RV $Y$ and $Z$ with i/o images $Iy(n,m)$, $Iz(n,m)$, also replace CDFs of $Y$, $Z$ with that of "computed from histograms". I denote it as $H_y(i)$ and $H_z(i)$ where $H_y(i)$ denotes the cumulative histogram obtained from image at hand, and $H_z(i)$ computed from the desired (to be matched) histogram. Under suitable conditions and correct scalings to ensure intensity levels to remain within valid limits, therefore we can make the definiton as $I_z = T(I_y)$ or more explicitly $I_z(n,m) = H_z^{-1}(H_y(I_y(n,m)))$

Your only problem is about the inverse cumulative histogram of $Z$ (or $I_z$) and unfortunately I don't know a method which may guarantee perfect inversion. As it is a CDF from a discrete histogram there might be empty bins and hence some indices of inverse $H_z(i)$ might have no mappings at all. It is up to you how to best overcome this problem.

However, as I understood from your code, you add pixels to empty/weak bins so that $H_z(i)$ becomes invertable...

Below I will put some matlab code "as is" : There are many methods to match histograms and I don't think below is the most accurate or the efficient one.

Iy = imread('cameraman.tif');  % Intesity between [0-255]
hy = imhist(Iy);         % compute image histogram hy(i)
figure,stem(hy);
title('original image histogram');

S = size(Iy);
Imax = 255;             % 8 bit
K = Imax / (S(1)*S(2)); % Scale factor
Hy = zeros(1,256);      % compute CDF Hy(i) of input image
Hy(1)=hy(1);
for i=2:256
Hy(i) = Hy(i-1) + hy(i);
end
Hy = K*Hy;
figure,stem(Hy);
title('CDF of original input image');

% Get the desired image and histograms:
figure,imshow(I_des);
S = size(I_des);

h_des = imhist(I_des);
figure,stem(h_des);
title('desired image histogram');

h_des(h_des < (S(1)*S(2)/(Imax*8))) = h_des(h_des < (S(1)*S(2)/(Imax*8))) + 32; % HERE  paddings
Hz = zeros(1,256);  % This is desired CDF
Hz(1)= h_des(1);
for i=2:256
Hz(i)=Hz(i-1)+h_des(i);
end
K_des = Imax / (Hz(256));
Hz = K_des*Hz;
figure,stem(Hz);
title('Desired CDF Hz');

% MATCH is performed
Imatch = (Iy);     % Processed image from matching
for i=0:255       % for each intensity level
ind = (Iy==i);   % find those index with an intensity value of "i"

j = Hy(i+1);

looping=1;  %  implement you own method
while looping
k = find( (Hz > (j-1)) & (Hz < (j+1)));
if length(k)>0   % for practical matters...
looping =0;
else
j = j + 1;    % when there is no match, update a little
end
end

Imatch(ind) = k(1); % Adjust the original intensity now...
end

figure, imshow(Imatch);
title('Image Matched to desired histogram');

h_match = imhist(Imatch);
figure, stem(h_match);
title(' Matched Histogram, not a good match ???');

• Thanks for the sound mathematical derivation! Helped me a lot to understand the problem conceptually. As for calculating the inverse histogram, I observed that it didn't cause much problem since there were not many empty bins in the histogram. – chronosynclastic Jan 16 '15 at 11:33
• I also managed to re-write the Matlab code you have given in Python, the only problem now is that it takes a very long time for a single image, since my algorithm requires to do this operation for every single row of an image separately, thus, when take into account the for-loop for every intensity level, the number of iterations turns out to be: # (rows in the image) * 256. I will try to find a more efficient way to implement that though. – chronosynclastic Jan 16 '15 at 11:38
• Glad to see if it helped. For the efficiency problem, it takes less than a second with an intel atom processor in my testing. As you said, it matches histogram of the entire image, instead of per row histograms. Matlab has very efficient vectorization (such as the line "ind = (I==i)" ) I am not sure if this can be as efficiently done in python (may be more efficient?). I could only suggest vectorization instead of explicit loop calculations. Also isolate the bottleneck, and lets see if we can change the algorithm to overcome that. – Fat32 Jan 16 '15 at 21:01