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     
img = cv2.imread('images/flicker1.jpg',0)
rows,cols = img.shape
# Number of pixels to pad each bin of a row histogram
N = 30
########################################################################
cdf_hist_Padded = np.zeros((rows,256))
cdf_hist_noPad = np.zeros((rows,256))
for i in range(0,rows):    
    # Read one row
    img_row = img[i,]
    # Calculate the row histogram (without padding)
    hist_row_noPad = cv2.calcHist([img_row],[0],None,[256],[0,256])
    # Calculate the cumulative row histogram (without padding)
    cdf_hist_row_noPad = hist_row_noPad.cumsum()
    cdf_hist_noPad[i,:] = cdf_hist_row_noPad
    # Copy the row to prepare for padding
    hist_row_Padded = np.copy(hist_row_noPad)
    # Apply uniform padding to all the bins less than T
    hist_row_Padded[0:T] = hist_row_Padded[0:T] + N
    # Calculate the cumulative row histogram of the padded row
    cdf_hist_row_Padded = hist_row_Padded.cumsum()
    # Accumulate the cumulative histograms of padded rows in a matrix
    cdf_hist_Padded[i,:] = cdf_hist_row_Padded

# 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_cdf_hist_noPad = cdf_hist_noPad/cdf_hist_noPad.max()
norm_cdf_hist_Padded = cdf_hist_Padded/cdf_hist_Padded.max()
norm_Gauss_cdf = Gauss_cdf/Gauss_cdf.max()

# Take the first original and padded+smoothed row cumulative histogram 
H_zero = norm_cdf_hist_noPad[0,:]
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 &lt; 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&lt; y) = P(G(X)&lt;y) = P(X&lt; 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:
I_des = imread('tire.tif');
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 &lt; (S(1)*S(2)/(Imax*8))) = h_des(h_des &lt; (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 &gt; (j-1)) &amp; (Hz &lt; (j+1)));
    if length(k)&gt;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 ???');