基於小波變換的多聚焦映像融合演算法——Python實現

# 畢設要做的東西，之前在CSDN上查資料，小波變換這部分全是MATLAB，而且沒幾個寫的靠譜的。Python的pywavelet庫資料太少了，我就把自己做的東西傳上來吧。

# 什麼你問我為什麼不用MATLAB。嗯嗯，我討厭這個軟體，就這樣。

# -*- coding: utf-8 -*-

import pywtimport numpy as npimport matplotlib.pyplot as pltfrom PIL import Imageimport cv2'''論文中的兩種方案：1，對低頻分量的所有像素點計算其局部方差，每張圖所有點的方差加起來除以兩張圖所有點加起來，得到兩張圖的權重融合映像每個像素點的值為兩張圖對應像素點的值加權平均，這個權就是上面算出來的權值。2，對每個高頻分量的像素點用canny運算元進行邊緣提取，再對邊緣映像的每一個像素點計算其局部方差，得到方差映像融合映像每個像素點的值為方差圖片中對應像素點取值較大的那張圖片的像素點。3，效果不是很好，在小波分塊的邊緣有明顯的灰階跳變(就是有些論文裡說的分塊效應)，but why?4，已找到原因，要求的不是某一個點對全圖的方差，而是在某點附近開個小視窗求視窗的局部方差5，現用局部方差的方法對多聚焦映像效果非常完美'''def imgOpen(path):    img=Image.open(path).convert('L')    imgArray=np.array(img)    return imgArray# 對於低頻分量，計算兩圖的權重比def varianceWeight(img1,img2):    mean1,var1=cv2.meanStdDev(img1)    mean2,var2=cv2.meanStdDev(img2)    weight1=var1/(var1+var2)    weight2=var2/(var1+var2)    return weight1,weight2# 實測這個函數效果非常好。。。def getVarianceImg(array):    row,col=array.shape    varImg=np.zeros((row,col))    for i in xrange(row):        for j in xrange(col):            up=i-5 if i-5>0 else 0            down=i+5 if i+5<row else row            left=j-5 if j-5>0 else 0            right=j+5 if j+5<col else col            window=array[up:down,left:right]            mean,var=cv2.meanStdDev(window)            varImg[i,j]=var    return varImg# 不會寫canny，暫時先用Sobel運算元代替def calcGradient(img):    xDiff=cv2.Sobel(img,cv2.CV_16S,1,0)    yDiff=cv2.Sobel(img,cv2.CV_16S,0,1)    stdXdiff=cv2.convertScaleAbs(xDiff)    stdYdiff=cv2.convertScaleAbs(yDiff)    gradient=np.sqrt(stdXdiff**2+stdYdiff**2)    return gradientdef testWave(img1,img2):    transf1=pywt.wavedec2(img1,'haar',level=4)    transf2=pywt.wavedec2(img2,'haar',level=4)    assert len(transf1)==len(transf2)    recWave=[]    for k in range(len(transf1)):        # 處理低頻分量        if k==0:            loWeight1,loWeight2 = varianceWeight(transf1[0],transf2[0])            lowFreq = np.zeros(transf2[0].shape)            row,col = transf1[0].shape            for i in range(row):                for j in range(col):                    lowFreq[i,j] = loWeight1*transf1[0][i,j] + loWeight2*transf2[0][i,j]            recWave.append(lowFreq)            continue        # 處理高頻分量        cvtArray=[]        for array1,array2 in zip(transf1[k],transf2[k]):            tmp_row,tmp_col = array1.shape            highFreq = np.zeros((tmp_row,tmp_col))            var1=getVarianceImg(array1);var2=getVarianceImg(array2)            for i in range(tmp_row):                for j in range(tmp_col):                    highFreq[i,j]=array1[i,j] if var1[i,j]>var2[i,j] else array2[i,j]            cvtArray.append(highFreq)        recWave.append(tuple(cvtArray))    return pywt.waverec2(recWave,'haar')def testPlot(org1,org2,img):    plt.subplot(131)    plt.imshow(org1,cmap='gray')    plt.axis('off')    plt.subplot(132)    plt.imshow(org2,cmap='gray')    plt.axis('off')    plt.subplot(133)    plt.imshow(img,cmap='gray')    plt.axis('off')    plt.show()if __name__=='__main__':    img1=imgOpen('F:\\Python\\try\\BasicImageOperation\\pepsia.jpg')    img2=imgOpen('F:\\Python\\try\\BasicImageOperation\\pepsib.jpg')    rec=testWave(img1,img2)    testPlot(img1,img2,rec)