【图像增强】梯度加权直方图均衡的对比度增强算法

梯度加权直方图均衡和传统直方图的比较

传统直方图均衡算法是一种用于图像增强的技术，其主要思想是通过对图像像素值的分布进行变换，将原始图像中像素值的分布变换为更均匀的分布，从而使得图像的视觉效果更好。该算法易于实现，并且可以显著提高对比度，使得图像中不同区域的颜色更加鲜明，更加易于观察和分析，而且通用性非常强。

但是该算法对噪声敏感，当图像中存在噪声时，会使得处理后的图像更加模糊，降低图像质量。此外，当图像中存在高亮或低亮区域时，直方图均衡算法容易使得这些区域过度饱和或过度降低亮度，从而降低图像的质量。直方图均衡算法会改变图像的像素值分布，有可能会引入新的问题，例如出现明显的条纹或噪点。当图像中存在大量相同像素的点且没有梯度变化时，很容易造成图像突变。

因此我们提出一种新的直方图均衡算法，需要考虑的是我们在灰度级拉伸的时候需要尽可能考虑到把梯度变化大的像素附近作大拉伸，梯度变化小的亮度值附近作小拉伸，我们对该算法起名为梯度加权直方图均衡算法

先看效果：

上面3幅图像中左图为原图，可以看到该x光图很暗而且背景有大片纯色区域，如果采用传统直方图均衡（中图）就会把大片纯色区域的亮度值抬高导致图像失真。而如果考虑了梯度信息（右图），则均衡后就能充分考虑到图像中哪里有突变，就能自适应地拉伸。

上面三幅直方图分别对应原三幅图像，可以看到考虑了梯度信息后，对于梯度小的区域几乎无拉伸（背景），对于含有丰富信息的手部关节则平坦拉伸，有很好的增强效果。

算法详解

对于灰度图片构成的矩阵 A = { a i j } m × n \mathbf{A} = \{ a_{ij} \}_{m \times n} A={aij​}m×n​ ，首先对其进行高斯模糊，然后利用卷积核（Sobel算子）
$${\mathbf{C}}_x=\left(\begin{array}{ccc}1& 0& -1\\ 2& 0& -2\\ 1& 0& -1\end{array}\right),{\mathbf{C}}_y=\left(\begin{array}{ccc}1& 2& 1\\ 0& 0& 0\\ -1& -2& -1\end{array}\right),$$

对图像卷积得到 G x = { g i j 1 } , G y = { g i j 2 } \mathbf{G}_x = \{g^1_{ij}\} , \mathbf{G}_y = \{g^2_{ij}\} Gx​={gij1​},Gy​={gij2​}

该操作实际上获得了图像中每个点的水平梯度和竖直梯度，将二者合并就是带有方向的梯度。这里我们只取其模长。
合并得到梯度阵 G = { g i j } \mathbf{G} = \{g_{ij}\} G={gij​}，其中
$$g_{ij}=\sqrt{(g_{ij}^1{)}^2+(g_{ij}^2{)}^2}$$

构造灰度级序列
$$T(n)=\sum _{a_{ij}=n}{g}_{ij}$$

该序列即表示了每个像素值都对应了哪些梯度，对其加权求和，权重就是梯度。比较传统直方图均衡的算法，其灰度级序列仅仅是 $T^{\prime }(n)={\sum }_{a_{ij}=n}1$ ,这并未考虑到梯度的影响，因此会造成对纯色背景与含有丰富细节的部分作等同的拉伸处理，导致信息丢失和对比度增强出现问题。

当然可以对梯度等级作离散分段处理。我们可以考虑使用边缘检测对超过可见范围的跳变像素值增强，其余的不变，这样对背景噪声的抑制效果更强。
设通过常规边缘检测算法（如Canny）得到边缘检测阵 E = { e i j } \mathbf{E} = \{ e_{ij} \} E={eij​}，边缘点为 $1$ 否则为 $0$ .
将梯度等级分为 $p$ 份（例如 $5$）：设等级分布阵梯度阵 L = { l i j } \mathbf{L} = \{ l_{ij} \} L={lij​}，其中权值序列：
$$l_{ij}=\text{ReLU}\left(\lceil \frac{p{g}_{ij}}{\max (\mathbf{G})}-1\rceil \right)$$
然后得到每个梯度等级的权重：
$${\gamma }_k=\frac{{\sum }_{l_{ij}=k}{e}_{ij}}{{\sum }_{l_{ij}}{e}_{ij}},0\le k\le p-1$$
设每个区域对应的权值阵为 W = { w i j } W = \{ w_{ij} \} W={wij​}，满足
$$w_{ij}={\gamma }_{l_{ij}}$$
因此构造灰度级序列也可以是如下形式：
$$T(n)=\sum _{a_{ij}=n}{w}_{ij}{e}_{ij}$$

归一化灰度级变换序列
$$S(n)=\frac{255}{{\sum }_{k=0}^{255}T(k)}\cdot \sum _{k=0}^{n}T(k)$$

该对灰度级序列的积分操作实际上就是对源图像的拉伸函数，将原图像像素经过上面公式变换后就得到了新图像。为了保证变换后的图像范围仍旧在 $0$ 到 $255$ 之间，所以加了一个拉伸因子 $\frac{255}{{\sum }_{k=0}^{255}T(k)}$ 用以保证映射范围。

均衡化得到的新图像 B = { b i j } m × n B = \{b_{ij}\}_{m \times n} B={bij​}m×n​ 满足
$$b_{ij}=\lfloor S(a_{ij})+0.5\rfloor$$

这一步 $+0.5$ 然后向下取整的操作是达到四舍五入的目的。

代码实现

给出python代码示例，实际上利用python的opencv库可以很容易实现

仅考虑梯度加权

import numpy as np
import cv2
def gradient_equalization(image):
    image = cv2.GaussianBlur(image, (5, 5), 0, borderType=cv2.BORDER_REPLICATE)
    sobel_x = cv2.Sobel(image, cv2.CV_64F, 1, 0, ksize=3)
    sobel_y = cv2.Sobel(image, cv2.CV_64F, 0, 1, ksize=3)
    gradients = np.sqrt(sobel_x**2 + sobel_y**2)
    GrayTrans = np.zeros(256)
    for i in range(256):
        GrayTrans[i] = GrayTrans[i - 1] + np.sum(gradients[image == i])
        
    GrayTrans = np.uint8(GrayTrans * 255 / GrayTrans[255] + 0.5)
    return GrayTrans[image]

考虑边缘检测的分级加权

def weighted_equalization(image, parts=5, threshold1=50, threshold2=100):
    sobel_x = cv2.Sobel(image, cv2.CV_64F, 1, 0, ksize=3)
    sobel_y = cv2.Sobel(image, cv2.CV_64F, 0, 1, ksize=3)
    gradients = np.sqrt(sobel_x**2 + sobel_y**2)
    edges = cv2.Canny(image, threshold1, threshold2)
    level = np.uint8(np.maximum(0, np.ceil(gradients * parts / np.max(gradients))-1))
    weight = np.zeros(parts)
    for i in range(parts):
        weight[i] = np.sum(edges * (level == i))

    weight /= np.sum(edges)
    GrayTrans = np.zeros(256)
    for i in range(256):
        GrayTrans[i] = GrayTrans[i - 1] + np.sum(weight[level[image == i]] * edges[image == i])
    
    GrayTrans = np.uint8(GrayTrans * 255 / GrayTrans[255] + 0.5)
    return GrayTrans[image]

毕竟python的循环是真的慢！！！！（无奈）对于大尺寸图象运行时间越来越长…可以对图像下采样之后得到灰度变换序列并用该序列去映射源图像，这样基本没有什么差别。

效果展示和结语

再展示一个个例子
左中右分别为：原图、传统直方图均衡、基于梯度直方图均衡，可以看到传统直方图均衡把天空噪声给放大了很多。

后面会更新基于累计边缘点可视范围的多级量化算法，该算法也是考虑到了最大程度展示原图像的细节，敬请期待。