二維影像haar小波變換的分解與重構

二維離散小波的理論推導和一維小波類似，但是以其尺度函式生成的尺度函式集作為標準正交基的尺度空間V_i的正交補空間W_i不能直接得到，而是可以證明，正交補空間W_i是由三個子空間的直和組成，對應的三個子空間可以由作為正交基的尺度函式、小波函式張成，

二維離散小波變換對影像的分解可以看做如下圖所示的濾波程序，即首先進行行濾波，沿著列方向進行，然后下采樣，然后對上一步得到的結果進行列濾波，沿著行方向，然后下采樣，做完所有的列濾波后，獲得4個不同的頻帶，一個近似分量、三個細節分量（水平、垂直、對角線），將所有的結果組合為一張圖，

若對所得的近似分量繼續進行這樣的濾波程序，既可以得到如下圖所示的塔式分解：

注意這里有個坑，python cv2對于float型別會乘以255，導致近似分量會顯示為白圖，需要使用pyplot進行顯示，重構出來的影像如下圖所示：

haar代碼如下：

import numpy as np
import matplotlib.image as mpimg # mpimg 用于讀取圖片
import matplotlib.pyplot as plt

L_D = np.array([0.7071,   0.7071])
H_D = np.array([-0.7071,  0.7071])
L_R = np.array([0.7071,   0.7071])
H_R = np.array([0.7071,  -0.7071])

QP = 20
threshold = 5
#這里的影像是整張影像，height、width是整張影像大小
def haar_decomposition(img, height, width, depth):
    img_h = np.zeros((height , width ))
    img_l = np.zeros((height , width ))

    img_l_l = np.zeros((height , width)) #近似系數，左上
    img_l_h = np.zeros((height , width)) #細節系數，垂直方向，左下
    img_h_l = np.zeros((height , width)) #細節系數，水平方向，右上
    img_h_h = np.zeros((height , width)) #細節系數，對角線，右下

    img_new = np.zeros((height , width)) 

    depth_count = 1;
    while depth_count <= depth:

        if depth_count != 1:
            width = width // 2
            height = height // 2
        
        #行濾波
        #行間濾波,然后進行下采樣
        for i in range(0, height // 2, 1):
            for j in range(0, width, 1):    
                i_endflag = 1 if i + 1 < height // 2 else -1
                
                img_l[i][j] = ((L_D[0] * img[i * 2][j] + L_D[1] * img[i * 2 + i_endflag][j]) // QP) * QP
                img_h[i][j] = ((H_D[0] * img[i * 2][j] + H_D[1] * img[i * 2 + i_endflag][j]) // QP) * QP
               
        n_w = width // 2
        n_h = height // 2
        #列濾波，列間濾波，沿著行
        for i in range(0, height // 2, 1):
            for j in range(0, width // 2, 1):
                j_endflag = 1 if j + 1 < width // 2 else -1
                img_l_l[i][j]             = ((L_D[0] * img_l[i][j * 2] + L_D[1] * img_l[i][j * 2 + j_endflag]) // QP) * QP
                img_l_h[n_h + i][j]       = ((H_D[0] * img_l[i][j * 2] + H_D[1] * img_l[i][j * 2 + j_endflag]) // QP) * QP
                img_h_l[i][n_w + j]       = ((L_D[0] * img_h[i][j * 2] + L_D[1] * img_h[i][j * 2 + j_endflag]) // QP) * QP
                img_h_h[n_h + i][n_w + j] = ((H_D[0] * img_h[i][j * 2] + H_D[1] * img_h[i][j * 2 + j_endflag]) // QP) * QP
    
        if depth_count == depth:
            img_new = img_new + img_l_l + img_l_h + img_h_l + img_h_h
        else:
            img_new = img_new + img_l_h + img_h_l + img_h_h
        
        img = img_l_l.copy()

        img_l[:][:] = 0
        img_h[:][:] = 0
        img_l_l[:][:] = 0    
        img_l_h[:][:] = 0
        img_h_l[:][:] = 0
        img_h_h[:][:] = 0
        
        depth_count += 1 
   
    img_new = np.round(img_new)
    
    return img_new
    
#這里的影像是整張影像，height、width是原始影像大小
def haar_reconstruction(img, height, width, trans_num):
    img_h = np.zeros((height , width ))
    img_l = np.zeros((height , width ))

    img_l_l = np.zeros((height , width))
    img_l_h = np.zeros((height , width))
    img_h_h = np.zeros((height , width))
    img_h_l = np.zeros((height , width))

    img_new = np.zeros((height , width)) 

    depth_count = 1;
    
    min_width  = width  // ( 2 ** trans_num )
    min_height = height // ( 2 ** trans_num )
    
    cur_width =  min_width
    cur_height = min_height
    
    while depth_count <= trans_num:

        if depth_count != 1:
            cur_width  = cur_width * 2
            cur_height = cur_height * 2
        
        #列上采樣
        for i in range(cur_height):
            for j in range(cur_width):
                img_l_l[i][2*j] = img[i][j]
                img_l_h[i][2*j] = img[cur_height + i][j]
                img_h_l[i][2*j] = img[i][cur_width + j]
                img_h_h[i][2*j] = img[cur_height + i][cur_width + j]
        
        #列濾波,列間濾波，沿著行
        for i in range(cur_height):
            for j in range(2 * cur_width):
                j_endflag = 1 if j + 1 < 2 * cur_width else -1
                
                img_l_l[i][j] = L_R[0] * img_l_l[i][j] + L_R[1] * img_l_l[i][j + j_endflag]
                img_l_h[i][j] = H_R[0] * img_l_h[i][j] + H_R[1] * img_l_h[i][j + j_endflag]
                img_h_l[i][j] = L_R[0] * img_h_l[i][j] + L_R[1] * img_h_l[i][j + j_endflag]
                img_h_h[i][j] = H_R[0] * img_h_h[i][j] + H_R[1] * img_h_h[i][j + j_endflag]
        
        for i in range(cur_height):
            for j in range(2 * cur_width):
                img_l[i][j] = img_l_l[i][j] + img_l_h[i][j]
                img_h[i][j] = img_h_l[i][j] + img_h_h[i][j]
        
                
        #行上采樣
        for i in range(cur_height - 1, -1, -1):
            for j in range(2 * cur_width):
                img_l[2 * i][j] = img_l[i][j]
                img_h[2 * i][j] = img_h[i][j]
                if i != 0:
                    img_h[i][j]     = 0
                    img_l[i][j]     = 0
                    
        #行濾波
        for j in range(2 * cur_width):
            for i in range(2 * cur_height):
                i_endflag = 1 if i + 1 < 2 * cur_height else -1
                img_l[i][j] = L_R[0] * img_l[i][j] + L_R[1] * img_l[i + i_endflag][j]
                img_h[i][j] = H_R[0] * img_h[i][j] + H_R[1] * img_h[i + i_endflag][j]

        img_new = img_l + img_h    
        img_new = np.round(img_new)
        img_new = np.clip(img_new, 0, 255)    
        
        for i in range(2 * cur_height):
            for j in range(2 * cur_width):
               img[i][j] = img_new[i][j]
        
        img_l[:][:] = 0
        img_h[:][:] = 0
        img_l_l[:][:] = 0    
        img_l_h[:][:] = 0
        img_h_l[:][:] = 0
        img_h_h[:][:] = 0
        
        depth_count += 1      
    
    return img_new

def calc_psnr(img_reco, img_org):
    data1 = np.round(img_reco)
    data2 = np.round(img_org)
    data1 = np.clip(data1, 0, 255)
    data2 = np.clip(data2, 0, 255)
    
    mse = np.mean((data1 - data2) * (data1-data2))
    # print(mse)
    return 10 * np.log10(255*255.0/mse)

# 統計零值占比
def calc_percent(img):
    count = 0
    height = img.shape[0]
    width = img.shape[1]
    for i in range(height):
        for j in range( width):
            if img[i, j] == 0:
                count += 1

    percent = count / 256 / 256
    return percent


def main():

    img = mpimg.imread("LENA.BMP")
    
    img = img.copy()
    trans_num = 1
    #img_st, 存盤小波變換后的影像
    img_st = haar_decomposition(img, img.shape[0], img.shape[1], trans_num)
    plt.xticks([])
    plt.yticks([])
    plt.imshow(img_st, cmap='gray')
    plt.show()

    per = calc_percent(img_st)
    
    img_re = haar_reconstruction(img_st, img.shape[0], img.shape[1], trans_num)
    mse = calc_psnr(img, img_re)

    plt.imshow(img_re, cmap='gray')
    plt.show()
    
    print(mse)
    print(per)
    
if __name__ == "__main__":
    main()