How to do histogram matching with normal distribution as reference?

Question

I have an image and I want to do histogram matching with normal distribution as reference. I've seen this link: Histogram matching with OpenCV, scikit-image, and Python. According to this link, the reference must be an image. But I don't have any images with normal distribution. How can I do this?

Answer 1

In general, we can use the answer from the following post, and replace the reference CDF with (normalized) cumulative sum of normally distributed data.

For computing the CDF of an image, we are collecting histogram, compute the cumulative sum, and normalize the result.

We may skip the image, and start computing the the CDF from the histogram.
Since we want to match with normal distribution, we may use normal distribution data as "histogram".

• Generate gaussian with normal distribution.

   # https://www.askpython.com/python/normal-distribution
   mean = 0.5
   sigma = 0.3  # Define STD to be 0.3 out of range [0, 1]
   x = np.linspace(0, 1, 256)
   prob_density = (np.pi*sigma) * np.exp(-0.5*((x-mean)/sigma)**2)
  

• Compute reference CDF from the normal distribution data:

   cum_t = np.cumsum(prob_density)
   c_t = cum_t / cum_t.max()  # Convert cumulative sum to CDF - divide by the maximum value.
  

• Execute cdf and hist_matching from the answer:

   c = cdf(img)
   matched_im = hist_matching(c, c_t, img)
  

---

Code sample:

from skimage.exposure import cumulative_distribution
from skimage import io
import matplotlib.pylab as plt
import numpy as np
from scipy.stats import norm

#https://stackoverflow.com/a/51702334/4926757
def cdf(im):
    '''
    computes the CDF of an image im as 2D numpy ndarray
    '''
    c, b = cumulative_distribution(im) 
    # pad the beginning and ending pixels and their CDF values
    c = np.insert(c, 0, [0]*b[0])
    c = np.append(c, [1]*(255-b[-1]))
    return c

def hist_matching(c, c_t, im):
    '''
    c: CDF of input image computed with the function cdf()
    c_t: CDF of template image computed with the function cdf()
    im: input image as 2D numpy ndarray
    returns the modified pixel values
    ''' 
    pixels = np.arange(256)
    # find closest pixel-matches corresponding to the CDF of the input image, given the value of the CDF H of   
    # the template image at the corresponding pixels, s.t. c_t = H(pixels) <=> pixels = H-1(c_t)
    new_pixels = np.interp(c, c_t, pixels)
    im = (np.reshape(new_pixels[im.ravel()], im.shape)).astype(np.uint8)
    return im

img = io.imread('Lenna_darken.png', as_gray=True)  # Read input image

# Generate gaussian with normal distribution.
# https://www.askpython.com/python/normal-distribution
mean = 0.5
sigma = 0.3  # Define STD
x = np.linspace(0, 1, 256)
prob_density = (np.pi*sigma) * np.exp(-0.5*((x-mean)/sigma)**2)

# Use cumulative sum of "prob_density" as replacement to CDF of the "template image".
cum_t = np.cumsum(prob_density)
c_t = cum_t / cum_t.max()  # Convert cumulative sum to CDF - divide by the maximum value.

c = cdf(img)
matched_im = hist_matching(c, c_t, img)

# Show image and plots for testing
plt.figure()
plt.imshow(matched_im, cmap='gray')
plt.figure()
plt.plot(x, prob_density)
plt.figure()
plt.plot(x, cum_t)
plt.figure()
plt.plot(x, c)
plt.figure()
plt.hist(matched_im.ravel(), bins=30)
plt.show()

io.imsave('matched_im.png', matched_im)

---

Input image:

Output image:

Normal distribution data:

Histogram of the output:

Note: result is not very accurate, due to quantization issues (working with 256 discrete values).