Histogam matching. This should be part of skimage

hist_match.py

def hist_match(im, ref_im):
    if np.squeeze(im).ndim == 2:
        im = np.expand_dims(im, axis=2)
        ref_im = np.expand_dims(ref_im, axis=2)
    out_im = np.empty_like(im)
    
    for chan_num in range(im.shape[2]):
        chan, ref_chan = im[...,chan_num], ref_im[...,chan_num]
        
        cdf, bin_centers = cumulative_distribution(chan)
        ref_cdf, ref_bin_centers = cumulative_distribution(ref_chan)
        chan_percentiles = np.interp(chan.flat, bin_centers, cdf)
        flat_out = np.interp(chan_percentiles, ref_cdf, ref_bin_centers)
        out = flat_out.reshape(chan.shape)
        out_im[..., chan_num] = out
        
    return np.squeeze(out_im)