Histogram matching (specification) in Python
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 ???');