数字图像处理02：直方图均衡化imhist函数的python实现

版权声明：本文为博主原创文章，未经博主允许不得转载。

---

数字图像处理02：直方图均衡化imhist函数的python实现

1、直方图均衡

直方图均衡就是将灰度比较集中一定范围的图像通过灰度映射/变换，使其灰度值平均分布在灰度级的范围内，以得到图像的更多细节，从而达到对比度增强的效果。

2、直方图均衡化的原理公式

1. 连续情况

其中Pr( r )和Ps(s)分别是原始图像灰度级r和变换后图像的灰度级s的概率密度函数PDF，L-1为图像的最大灰度级。上一个公式表示的是r的PDFPr( r )和s的PDFPs(s)的关系。给定一个输入，则Pr( r )已知，就可以求出s，继而求出Ps(s)。r和s的映射关系即为均衡化关系。由Pr( r )和Ps(s)可以得出输入图像的直方图和变换后的输出图像的直方图。


再经过以上两个公式的推导，可以得出结论Ps(s)=1/(L-1),可知该变换s=t( r )就是取均值。
2. 离散情况
对于图像来说，其实应该是离散情况。


其中MN为图像的像素总个数。以上两个公式是“1”中公式的离散情况，在实际对图像进行操作中也应该是离散情况下进行操作计算。

3、离散情况举例

注：该例子对应于冈萨雷斯.《数字图像处理(第三版)》中文翻译版，第76面例3.5。

上图中的数据经过‘2’中的公式计算得到经过变换的r对应的s值分别为：



对S0~S7经过四舍五入取整：

可以看原始图像数据像素集中在r0~r3 即灰度0~3，经过均衡化变换后，将0变换为1，1变换为3，2变换为5，3、4变换为6，5、6、7变换为7，这样就将灰度分散开。所谓均衡变换也就是灰度级的映射变换，即将灰度级密度大值的映射成其它原来灰度级密度小的值，以达到一种均衡化/平均的效果。从r到s的映射将集中的r分散到其他灰度，即将灰度值r变为其他灰度值，达到均衡化。

4、直方图均衡化计算步骤

其中第四步计算累积直方图是为了下一步(第五步)计算变换后对应的灰度值，而第五步中先加0.5再作INT(取整)操作即是四舍五入取整操作(此操作对于编程实现非常有用)。
对例3.5的具体步骤如下：

注：此上两张图来源于许录平.科学出版社. 《数字图像处理》. 配套教学PPT。

5、Matlab中直方图均衡化的imhist函数

h=imhist(f, b)，其中f为输入图像，h为直方图，b是用来形成直方图的“统计堆栈”的数目。“统计堆栈”仅仅是灰度的一小部分。例如我们处理一幅uint8类的图像且设b=2，然后灰度范围被分成两部分：0至127和128至255。所得的直方图将有两个值：h(1)等于图像在[0,127]间隔内的像素数，h(2)等于图像在[128,255]间隔内的像素数。通过下列表达式就可以得到归一化的直方图：p=imhist(f, b)/numel(f)，numel(f)函数可以给出数组f中元素的个数(即图像中的像素数)。

6、python代码实现

1. 显示图像标题所需的字体

from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']

2. 需要导入的库

import numpy as np
import matplotlib.pyplot as plt
from scipy import misc

3.画条形图/直方图或者可视化所需要的包和函数创建。
由于用matplotlib自带的绘制条形图/直方图的包imhist绘图时出现了问题，所以就选择了之前用过的数据可视化的绘图包pyecharts，pyecharts函数包的是一个可以根据数据绘制条形图/折线图/饼状图的python绘图包，具体用法请见主页：pyecharts 文档。

# 画条形图/直方图 或者可视化。导入pyecharts函数包
import os
from pyecharts import Bar,Pie,Line
# 创建直方图/条形图绘制函数
def get_charts(x, y, label, type):
    
    if type==1:
        c = Pie('饼状图')
    elif type==2:
        c = Bar('条形图')
    elif type==3:
        c = Line('折线图')
    c.add(label, x, y, is_more_utils=True)
    c.show_config()
    c.render()
    os.system(r"render.html")

其中x和y分别是要绘制图形的x轴和y轴的数据，label为图形的标注/标签，type为1、2、3分别表示 ‘饼状图’ 、‘条形图’、‘折线图’，操作实现非常简单，而且绘制出来的图形非常好看。

4.imhist函数实现

# imhist函数构建，img为原始图像，L为其灰度级
def imhist(img, L):
    
    f = misc.imread(img)
    plt.figure(1)
    plt.imshow(f, cmap='gray')
    plt.axis('off')
    plt.title('原始图像')
    plt.show()
    
    
    w, h = f.shape
    f1 = np.zeros([w, h])
    pallel = dict((i,0) for i in range(L))
    # 计算原始图像灰度级像素个数ni
    for x in range(0,w):
        for y in range(0,h):
            i = f[x, y]
            if i in pallel:
                pallel[i] = pallel[i] + 1
            else:
                pallel[i] = 1
    
    key = [i for i in range(L)]     # key为灰度级，即0~255
    value1 = [pallel[i] for i in range(L)]      # value1为各灰度级的像素个数ni
    n = w * h
    
    #计算原始图像的直方图，value2为原始图像的直方图P(i) = ni/n = ni/(w*h)
    value2 = [value1[i]/n for i in range(L)]        
    
    # 计算累积直方图，value3为累积直方图Pj
    value3 = [0 for i in range(L)]             
    value3[0] = value2[0]
    for j in range(1, L):
        value3[j] = value3[j-1] + value2[j]
    
    # 计算灰度值变换，value4为由r变换为s后的灰度值
    value4 = [0 for i in range(L)]          
    for j in range(L):
        value4[j] = int((L-1) * value3[j] + 0.5)
    
    # 计算变换后各灰度级的像素个数nj，value5即为nj
    value5 = [0 for j in range(L)]
    for x in range(0,w):
        for y in range(0,h):
            i = f[x, y]
            t = value4[i]
            f1[x, y] = t
            for k in range(L):
                if (t == k):
                    value5[k] = value5[k] + 1
#            j = f1[x, y]
#            if t in value5:
#                value5[t] = value5[t] + 1
#            else:
#                value5[t] = 1

    
    # 计算变换后图像的直方图，value6即为P(j)
    value6 = [value5[i]/n for i in range(L)]
    
    
    # 均衡化后的图像
    scipy.misc.imsave('huafen2_2.png', f1)
    plt.figure(2)
    plt.imshow(f1, cmap='gray') 
    plt.axis('off')
    plt.title('变换图像')
    plt.show()
    
    return f1, f, pallel, key, value1, value2, value3, value4, value5, value6

(1) 其中，img为输入图像，L为最大灰度级，如256。

(2)

    f = misc.imread(img)
    plt.figure(1)
    plt.imshow(f, cmap='gray')
    plt.axis('off')
    plt.title('原始图像')
    plt.show()

这段代码用来读入原始图像。

(3)

    pallel = dict((i,0) for i in range(L))
    # 计算原始图像灰度级像素个数ni
    for x in range(0,w):
        for y in range(0,h):
            i = f[x, y]
            if i in pallel:
                pallel[i] = pallel[i] + 1
            else:
                pallel[i] = 1
    
    key = [i for i in range(L)]     # key为灰度级，即0~255
    value1 = [pallel[i] for i in range(L)]      # value1为各灰度级的像素个数ni
    n = w * h

pallel表示由原图像的灰度值和该灰度值的个数作为键值对的字典，灰度值为键，个数为值。首先创建一个长度为L的空字典，然后遍历图像，通过两个判断语句if和else将图像灰度值的个数添加到字典pallel中，如果灰度值i在pallel中，则pallel[i]对应的i这个灰度值的个数加一，否则 灰度值i的个数等于一，依次遍历整个图像矩阵。最后将灰度级0~255存到列表key中，将各灰度级的个数存到value1中。

(4) 然后按照 ‘4、’ 中的8个步骤依次计算剩余的变量/数据列表。

#计算原始图像的直方图，value2为原始图像的直方图P(i) = ni/n = ni/(w*h)
    value2 = [value1[i]/n for i in range(L)]        
    
    # 计算累积直方图，value3为累积直方图Pj
    value3 = [0 for i in range(L)]             
    value3[0] = value2[0]
    for j in range(1, L):
        value3[j] = value3[j-1] + value2[j]
    
    # 计算灰度值变换，value4为由r变换为s后的灰度值
    value4 = [0 for i in range(L)]          
    for j in range(L):
        value4[j] = int((L-1) * value3[j] + 0.5)
    
    # 计算变换后各灰度级的像素个数nj，value5即为nj
    value5 = [0 for j in range(L)]
    for x in range(0,w):
        for y in range(0,h):
            i = f[x, y]
            t = value4[i]
            f1[x, y] = t
            for k in range(L):
                if (t == k):
                    value5[k] = value5[k] + 1
#            j = f1[x, y]
#            if t in value5:
#                value5[t] = value5[t] + 1
#            else:
#                value5[t] = 1

    
    # 计算变换后图像的直方图，value6即为P(j)
    value6 = [value5[i]/n for i in range(L)]

(5) 保存并显示均衡化变换后代图像，并返回各数据列表

    # 均衡化后的图像
    scipy.misc.imsave('huafen2_2.png', f1)
    plt.figure(2)
    plt.imshow(f1, cmap='gray') 
    plt.axis('off')
    plt.title('变换图像')
    plt.show()
        
    return f1, f, pallel, key, value1, value2, value3, value4, value5, value6

5. 测试imhist函数，并绘制条形图/直方图

# imhist函数测试        
f1, f,pallel, key, value1, value2, value3, value4, value5, value6 = imhist('huafen2.tif', 256)

# 绘制原始图像灰度级像素个数的条形图
label1 = "灰度级:像素个数"
type = 2     # type = 2表示条形图
get_charts(key, value1, label1, type)

# 绘制原始图像灰度级概率密度的条形图
label2 = "灰度级:概率密度"
type = 2     # type = 2表示条形图
get_charts(key, value2, label2, type)

# 绘制原始图像的累积直方图
label3 = "灰度级:累积概率密度"
type = 2     # type = 2表示条形图
get_charts(key, value3, label3, type)

# 绘制映射变换关系的条形图
label3 = "灰度级:灰度级"
type = 2     # type = 2表示条形图
get_charts(key, value4, label3, type)

# 绘制变换后图像灰度级像素个数的条形图
label3 = "灰度级:像素个数"
type = 2     # type = 2表示条形图
get_charts(key, value5, label3, type)

# 绘制变换后图像灰度级概率密度的条形图
label3 = "灰度级:概率密度"
type = 2     # type = 2表示条形图
get_charts(key, value6, label3, type)

输入图像f为 冈萨雷斯.《数字图像处理(第三版)》中文翻译版 第三章的例图 “huafen”:

均衡化变换后的图像f1为：

经过均衡化变换后，原始的暗色图象个变亮了，可以得到图像的更多细节。
(1) value1的条形图(原始图像各灰度值的个数统计图)如下：

在该条形图中可以看出原始图像的各像素值的个数，以及其主要分布在较低的灰度级，所以原始图像看上去比较暗。
(2) value2的直方图(原始图像各灰度级的概率密度图)如下：

(3) value3的条形图(原始图像各灰度级的原始图像的累积直方图)如下：

(4) value4的条形图(原始图像各灰度级的映射关系图)如下：

从这幅图可以看到原始图像灰度值到变换后图像灰度值的映射关系。
(5) value5的条形图(变换后图像各灰度值的个数统计图)如下：

从这幅图可以看出，均衡化变换后图像的灰度值均匀的分布在灰度空间，达到了均衡化的效果。

(6) value6的直方图(变换后图像各灰度级的概率密度图)如下：

6. 整体代码

# -*- coding: utf-8 -*-
"""
Created on Sun Jan 20 09:20:58 2019

@author: ChengGD
"""

# 图像显示窗口标题需要的中文字体
from pylab import mpl
mpl.rcParams['font.sans-serif'] = ['SimHei']


#from PIL import Image
import numpy as np
import matplotlib.pyplot as plt
from scipy import misc
import scipy 
#import cv2


# 画条形图/直方图 或者可视化。导入pyecharts函数包
import os
from pyecharts import Bar,Pie,Line
# 创建直方图/条形图绘制函数
def get_charts(x, y, label, type):
    
    if type==1:
        c = Pie('饼状图')
    elif type==2:
        c = Bar('条形图')
    elif type==3:
        c = Line('折线图')
    c.add(label, x, y, is_more_utils=True)
    c.show_config()
    c.render()
    os.system(r"render.html")


# imhist函数构建，img为原始图像，L为其灰度级
def imhist(img, L):
    
    f = misc.imread(img)
    plt.figure(1)
    plt.imshow(f, cmap='gray')
    plt.axis('off')
    plt.title('原始图像')
    plt.show()
    
    
    w, h = f.shape
    f1 = np.zeros([w, h])
    pallel = dict((i,0) for i in range(L))
    # 计算原始图像灰度级像素个数ni
    for x in range(0,w):
        for y in range(0,h):
            i = f[x, y]
            if i in pallel:
                pallel[i] = pallel[i] + 1
            else:
                pallel[i] = 1
    
    key = [i for i in range(L)]     # key为灰度级，即0~255
    value1 = [pallel[i] for i in range(L)]      # value1为各灰度级的像素个数ni
    n = w * h
    
    #计算原始图像的直方图，value2为原始图像的直方图P(i) = ni/n = ni/(w*h)
    value2 = [value1[i]/n for i in range(L)]        
    
    # 计算累积直方图，value3为累积直方图Pj
    value3 = [0 for i in range(L)]             
    value3[0] = value2[0]
    for j in range(1, L):
        value3[j] = value3[j-1] + value2[j]
    
    # 计算灰度值变换，value4为由r变换为s后的灰度值
    value4 = [0 for i in range(L)]          
    for j in range(L):
        value4[j] = int((L-1) * value3[j] + 0.5)
    
    # 计算变换后各灰度级的像素个数nj，value5即为nj
    value5 = [0 for j in range(L)]
    for x in range(0,w):
        for y in range(0,h):
            i = f[x, y]
            t = value4[i]
            f1[x, y] = t
            for k in range(L):
                if (t == k):
                    value5[k] = value5[k] + 1
#            j = f1[x, y]
#            if t in value5:
#                value5[t] = value5[t] + 1
#            else:
#                value5[t] = 1

    
    # 计算变换后图像的直方图，value6即为P(j)
    value6 = [value5[i]/n for i in range(L)]
    
    
    # 均衡化后的图像
    scipy.misc.imsave('huafen2_2.png', f1)
    plt.figure(2)
    plt.imshow(f1, cmap='gray') 
    plt.axis('off')
    plt.title('变换图像')
    plt.show()
    
    return f1, f, pallel, key, value1, value2, value3, value4, value5, value6
    
# imhist函数测试        
 f1, f, pallel, key, value1, value2, value3, value4, value5, value6 = imhist('huafen2.tif', 256)

# 绘制原始图像灰度级像素个数的条形图
label1 = "灰度级:像素个数"
type = 2     # type = 2表示条形图
get_charts(key, value1, label1, type)

# 绘制原始图像灰度级概率密度的条形图
label2 = "灰度级:概率密度"
type = 2     # type = 2表示条形图
get_charts(key, value2, label2, type)

# 绘制原始图像的累积直方图
label3 = "灰度级:累积概率密度"
type = 2     # type = 2表示条形图
get_charts(key, value3, label3, type)

# 绘制映射变换关系的条形图
label3 = "灰度级:灰度级"
type = 2     # type = 2表示条形图
get_charts(key, value4, label3, type)

# 绘制变换后图像灰度级像素个数的条形图
label3 = "灰度级:像素个数"
type = 2     # type = 2表示条形图
get_charts(key, value5, label3, type)

# 绘制变换后图像灰度级概率密度的条形图
label3 = "灰度级:概率密度"
type = 2     # type = 2表示条形图
get_charts(key, value6, label3, type)