直方图均衡化与规定化

1. 直方图均衡化

一幅图像的灰度级可以看成是区间[0, L-1]内的随机变量。基本描绘子是灰度级的概率密度函数(PDF)。
令r为输入的灰度变量，s为输出的灰度变量，p_r(r) 和 p_s(s) 分别表示灰度随机变量r和s的概率密度函数，需求灰度变量 r -> s 的映射关系。若满足可微条件，则输出灰度变量s的PDF可由下面公式得到：

\[p_{s}(s)=p_{r}(r)\left|\frac{\mathrm{d} r}{\mathrm{d} s}\right| \]

在图像处理中特别重要的变换函数有如下形式：

\[s=T(r)=(L-1) \int_{0}^{r} p_{r}(w) \mathrm{d} w \]

其中w是积分假变量。公式右边是输入灰度随机变量r的累积分布函数(CDF)。因为PDF总为正，一个函数的积分是该函数下方的面积，而积分值为1（PDF曲线下方的面积总是1），积分区间为[0, L-1]。

对于一个均衡化后的图像，其概率密度均匀，即 p_s(s) =1 / (L-1) 。
则变换后的灰度变量 s 与变换前的灰度变量 r 的关系可由上式求出。

对于离散值，一幅图像中灰度级r_k出现的概率近似为：

\[p_{r}\left(r_{k}\right)=\frac{n_{k}}{M N}, \quad k=0,1,2, \cdots, L-1 \]

对应上述变换离散形式为：

\[s_{k}=T\left(r_{k}\right)=(L-1) \sum_{j=0}^{k} p_{r}\left(r_{j}\right)=\frac{(L-1)}{M N} \sum_{j=0}^{k} n_{j}, \quad k=0,1,2, \cdots, L-1 \]

其中MN是图像中像素的总数，n_k是灰度为r_k的像素个数，L是图像中灰度级数量。

Python代码实现一张16bit图的均衡化：

• 1.读入一张16bit图像矩阵，可以使用 img = cv2.imread('origin.tif', cv2.IMREAD_UNCHANGED) 和 img_array = np.array(img)方式读取；
• 2.使用np.histogram()计算图像直方图；
• 3.计使用np.cumsum()计算累计直方图；
• 4.将每个灰度级出现的次数除以 所有灰度级出现的总次数（即图像的宽x高）；
• 5.计算 63336 x 63336 维的灰度变换关系矩阵，把 r 的每个灰度级 变换到 s 的每个灰度级；
• 6.将整幅图每个像素根据上面的灰度变换关系矩阵进行灰度变换；
• 7.输出灰度变换完后的图像矩阵，即完成直方图均衡化

def histEqualize(imgArray):
    """计算16bit图进行直方图均衡化后的结果
    
    Args:
        imgArray: 原图灰度矩阵
        
    Returns:
        equalimgArray：均衡化后的图片灰度矩阵
    
    """
                
    imgHist = np.histogram(imgArray, bins=range(0, 65537))[0]
    # 累计直方图    
    imgCumuHist = np.cumsum(imgHist)
    # 累计直方图概率
    imgCumuHistPro = imgCumuHist / imgCumuHist[-1]
    # 灰度变换矩阵
    grayConvert = np.ceil(np.dot(65535, imgCumuHistPro))
    w,h = imgArray.shape
    equalimgArray = np.zeros((w,h), dtype=np.uint16)
    for i in range(w):
        for j in range(h):
            equalimgArray[i][j] = grayConvert[imgArray[i][j]]
   
    return equalimgArray

2. 直方图规定化（匹配）

令r和z分别表示输入图像和输出图像（已直方图规定化处理）的灰度级。给定输入图像灰度级概率密度函数为p_r(r) ，而 p_z(z) 是我们希望输出图像所具有的指定概率密度函数， 因此需要求 r -> z 灰度变量映射关系。
我们先将r进行直方图均衡变换 T(r) 为中间变量s:

\[s=T(r)=(L-1) \int_{0}^{r} p_{r}(w) \mathrm{d} w \]

再将z进行直方图均衡变换G(z)，也映射到中间变量s:

\[G(z)=(L-1) \int_{0}^{z} p_{z}(t) \mathrm{d} t=s \]

则有三个变量的关系为：

\[z=G^{-1}[T(r)]=G^{-1}(s) \]

因此直方图规定化步骤如下：

• 1.计算图像的直方图p_r(r)，并进行直方图均衡变换。计算 r -> s 的映射关系，把s四舍五入为范围[0, L-1]内的整数；
• 2.对配准变量z进行直方图均衡化，计算变换后G(z)的所有值。把G(z)的值四舍五人为范围[0,L-1]内的整数；
• 3.对于灰度变量s的每一个值s_k ,k=0,1.2,...,L-1, 使用步骤2存储的G(z)值寻找相应的z值,以使G(z)最接近s_k，并存储这些 s -> z 的映射。

Python使用numpy和opencv实现上述步骤如下。
其中在计算 s->z 的映射关系时，将s的所有值转变为一个集合set，避免对于同一s值，重复寻找最接近同时z最小的G(z)值。
而把其它灰度级置零z_Gray = np.zeros(65536)，是因为[0, 65535]的灰度级中，没有在s中出现的值，是不需要考虑的，因为不会出现在 r -> s -> z 的映射关系中。
寻找s与G(z)最接近的值的索引z的方式，是先通过G(z)列表减去全为s值的同维度列表，再寻找最小值的索引。

总结就是，输入图像r映射到s，再把配准变量z映射到G(z)=s，通过离散变量查找G(z)=s映射到z的关系。最终通过 r - > s -> z 完成灰度值映射。

def histSpecify(imgArray, specifyHist):
    """对16bit图进行直方图规定化（直方图匹配）
    
    Args:
        imgArray: 原图灰度矩阵
        specifyHist: 需要匹配的直方图
        
    Returns:
        specifiedImgArray: 规定化后的图片灰度矩阵
        
    """     
    
    # 原图的直方图
    imgHist = np.histogram(imgArray, bins=range(0, 65537))[0]  
    # 累计直方图    
    imgCumuHist = np.cumsum(imgHist)
    # 累计直方图概率
    imgCumuHistPro = imgCumuHist / imgCumuHist[-1]
    # 计算 r->s 灰度转换关系
    s_Gray = np.ceil(np.dot(65535, imgCumuHistPro))
    
    
    # 匹配直方图的累计图    
    specifyCumuHist = np.cumsum(specifyHist)
    specifyCumuHistPro = specifyCumuHist / specifyCumuHist[-1]
    # 计算 z->G(z) 灰度转换关系
    Gz_Gray = np.ceil(np.dot(65535, specifyCumuHistPro))
    
    # 计算 s->z 灰度转换关系
    z_Gray = np.zeros(65536)
    set_s_Gray = set(s_Gray)
    for s in set_s_Gray:
        # 寻找s与G(z)最接近的值的索引z，并作映射 s->z
        diff = np.abs(Gz_Gray - np.full(65536, s))
        z = np.argwhere(diff == min(diff))[0][0]
        z_Gray[int(s)] = z 
    
    # 图像灰度级变换
    w, h = imgArray.shape        
    specifiedImgArray = np.zeros((w,h), dtype = np.uint16)
    for i in range(w):
        for j in range(h):
             # r -> s
             specifiedImgArray[i][j] = s_Gray[imgArray[i][j]]
             # s -> z
             specifiedImgArray[i][j] = z_Gray[specifiedImgArray[i][j]]
    
    return specifiedImgArray

3. 图片测试

if __name__=='__main__':

    ori_img = cv2.imread('origin.tif', cv2.IMREAD_UNCHANGED)
    ori_img_array = np.array(ori_img)
    ori_hist = np.histogram(ori_img_array, bins=range(0, 65537))[0]
    # 在IPython中 %varexp --plot ori_hist 查看原图的直方图
       
    spec_img = cv2.imread('processed.tif', cv2.IMREAD_UNCHANGED)
    spec_img_array = np.array(spec_img)
    spec_hist = np.histogram(spec_img_array, bins=range(0, 65537))[0]
    # 在IPython中 %varexp --plot spec_hist 查看原图的直方图
    
    
    # %% cv方式显示直方图均衡化后的图片    
    equal_img_array = histEqualize(ori_img_array)
    cv2.imwrite('histeq_img.png', equal_img_array)
    img_scaled = cv2.normalize(equal_img_array, dst=None, alpha=0, beta=65535, norm_type=cv2.NORM_MINMAX)
    cv2.imshow('equal_img_array', img_scaled)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
        
    # %% cv方式显示直方图规定化后的图片    
    specified_img_array = histSpecify(ori_img_array, spec_hist)
    cv2.imwrite('specified_img.png', specified_img_array)
    img_scaled = cv2.normalize(specified_img_array, dst=None, alpha=0, beta=65535, norm_type=cv2.NORM_MINMAX)
    cv2.imshow('specified_img_array', img_scaled)
    cv2.waitKey(0)
    cv2.destroyAllWindows()

OpenCV方式计算图片的直方图

cv2.calcHist(images, channels, mask, histSize, ranges[, hist[, accumulate ]]) ->hist

• images:输入的图像

• channels:选择图像的通道

• mask:掩膜，是一个大小和image一样的np数组，其中把需要处理的部分指定为1，不需要处理的部分指定为0，一般设置为None，表示处理整幅图像

• histSize:使用多少个bin，一般为256，16bit图为65535

• ranges:像素值的范围，一般为[0,255]表示0~255，16bit图为[0,65535]

ori_hist = cv2.calcHist([ori_img], [0], None, [65536], [0, 65535])

GDAL库读取图片

ori_img = cv2.imread('origin.tif', cv2.IMREAD_LOAD_GDAL)