直方圖均衡化與規定化

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)