第十五節、韋伯局部描述符(WLD，附源碼)

本文轉載自查看原文 2018-09-01 19:29 1782 OpenCV

紋理作為一種重要的視覺線索，是圖像中普遍存在而又難以描述的特征，圖像的紋理特征一般是指圖像上地物重復排列造成的灰度值有規則的分布。紋理特征的關鍵在於紋理特征的提取方法。目前，用於紋理特征提取的方法有很多，最具有代表性的是有基於二階概率密度的灰度共生矩陣、符合人眼視覺特性的小波變換、紋理譜法以及基於圖像結構基元的紋理元方法等。
為了更有效地描述圖像局部紋理特征，又先后提取了局部二值模式(Local Binary Pattern,LBP)和韋伯局部描述符(Weber Local Descriptor,WLD),這兩種方法都是基於鄰域像素點間的灰度變化特征來描述圖像紋理的。因兩者易於理解、便於計算且具有較好的局部特征描述能力而被廣泛地應用於紋理分類、目標檢測、人臉識別、合成孔徑雷達(Synthetic Aperture Radar,SAR) 圖像索引、指紋活性檢測、圖像偽造檢測等眾多領域中。

一 WLD原理

韋伯定律是反映心理量和物理量之間關系的定律，它表明能夠引起感覺差異的差別閾值與原始刺激的強度之比是一個常量，即：

$$\frac{ΔI}{I}=k$$

其中$k$是一個常量；$ΔI$表示差別閾值；$I$表示原始刺激的強度。由此可以推知，刺激的變化所引起的感覺差異不僅與刺激變化的大小有關，還與原始刺激的強度有關。局部圖像描述符WLB就是根據該定律提出的，它包含兩個算子：差分激勵算子和方向算子。WLD計算除邊緣像素點外的每個像素點的差分激勵和方向，並以其二維分布直方圖來聯合表征圖像的紋理特征。

1、差分激勵算子

差分激勵反映局部窗內灰度變化的強度信息。通過計算局部窗內鄰域像素點與中心像素點間的灰度差值和與中心像素點灰度值的比值$G_{ratio}(x_c)$，再利用反正切變換將分布在$[-P,+∞]$范圍內的$G_{ratio}(x_c)$映射到區間$(-\frac{\pi}{2},\frac{\pi}{2})$內得到，差分激勵$ξ(x_c)$的計算式為：

$$ξ(x_c)=arctan(G_{ratio}(x_c))=arctan(\sum\limits_{i=0}^{P-1}\frac{(x_i-x_c)}{x_c})$$

其中,$x_c$和$x_i,i=0,1,...,P-1$分別表示中心像素點和鄰域像素點的灰度值；$P$表示鄰域像素點個數。

2、方向算子

方向反映局部窗內灰度變化的空間分布信息。通過局部窗內水平方向與垂直方向上鄰域像素點的灰度差值比值的反正切變換來描述。其計算式為:

$$Φ(x_c)=arctan(\frac{D_V}{D_H})$$

其中,$D_H$和$D_V$分別表示水平方向上和垂直方向上中心像素點兩側的鄰域像素點間的灰度差異，若對於$3\times3$像素的局部窗口如下圖所示，則$D_H=x_7-x_3$，$D_V=x_5-x_1$。

為了能夠更加有效的區分局部窗口的灰度分布變換，進一步將方向由$Φ(x_c)\in(-\frac{\pi}{2},\frac{\pi}{2})$變換到了$Φ'(x_c)\in[0,2\pi]$，其變換公式為：

$$Φ'(x_c)=\begin{cases} Φ(x_c) \qquad D_H>0,D_V>0 \\ \pi+Φ(x_c) \qquad D_H<0,D_V>0 \\\pi+Φ(x_c) \qquad D_H<0,D_V<0 \\ 2\pi+Φ(x_c) \qquad D_H>0,D_V<0 \end{cases}$$

3、WLD直方圖

WLD采用均勻量化技術，將方向$Φ'(x_c)$均勻地量化為$T$個方向，將差分激勵均勻地划分為$M$個頻段，分別對應於圖像中的高頻、中頻和低頻變化，再將划分的每個頻段上將差分激勵均勻地量化為$S$格，形成一個$T\times{C}=T\times(M\times{S})$的二維直方圖，並通過編碼將其轉化為一維向量用於表示圖像的紋理特征。

二維直方圖如上圖所示，橫軸表示方向，縱軸表示差分激勵。每個小矩形表示在該方向下所在差分激勵區間像素的數量，數量不同，顏色不同。

下圖是WLD直方圖得到的過程：

1、在每個主方向上得到差分激勵子直方圖，得到$H(0)$至$H(T-1)$;

2、將$H(k)$分成$M$個子區間，即$l_m,m=0,...,M-1$,將$H(k)$中的${l_m}$對應放置在$t=k$和$m=i$處；

3、將$m=i$的一行子直方圖拼接成一個直方圖，即$H_i$；

4、將$H_i$組合成一個直方圖，即WLD直方圖。

然后可以把WLD直方圖用於后續分類。

4、缺點

我們來看下圖所示的3個局部灰度分布示例。從3個局部窗內的灰度分布看，(a)~(c)分別表示了高頻、中頻和低頻3種變換。按照WLD的計算方法，他們的$ΔI$都等於0，即差分激勵等於0，而且，在垂直方向上的灰度值差值也為0，表明方向也等於0.這意味着，WLD將無法區分這3個紋理模式。

存在上述問題的主要原因有:

WLD的差分激勵算子是利用各向同性的邊緣檢測濾波器———拉普拉斯算子統計局部窗內$P$個鄰域像素點與中心像素點間的灰度差值之和$ΔI$，導致了灰度變化的正負差值相互抵消,換言之,局部窗內的灰度變化信息沒有充分體現;
WLD 的方向算子僅表達了水平方向和垂直方向上鄰域像素點間灰度變化梯度的空間分布方位,不能充分反映局部窗內灰度變化的空間分布結構信息,難以體現紋理的內在變化特征;

5、改進

紋理特征是指與空間分布相關的圖像灰度等級的變化。這意味着灰度圖像的紋理既與各像素點間灰度變化的梯度幅值有關,也與其梯度的空間分布密切相關,兩者是有機的一體。針對WLD 存在的問題，提出一種基於正負梯度改進的 WLD(WLD-PNG)。

其核心思想是:

基於局部窗內鄰域像素點與中心像素點間灰度變化的正負梯度信息,分離計算正、負差分激勵以保留灰度等級的變化特征;

利用 LBP 模式提取正負梯度分布的結構信息,以反映灰度等級變化的空間分布特征。最后,采用均勻量化和編碼技術,將兩者有機聯合建立圖像的紋理特征。

三代碼實現

由於opencv沒有提供WLD算法的實現，但是從網上我們可以找到基於Weber local descriptors的人臉識別的代碼和基於OPENCV的WLD特征提取程序，我們作為參考，修改自己的代碼如下：

# -*- coding: utf-8 -*-
"""
Created on Sat Sep  1 19:08:10 2018

@author: zy
"""

'''

WLD
韋伯局部描述符

'''
import cv2
import numpy as np


    
class WLD(object):
    def __init__(self):
        #差分激勵由中心像素與周圍像素強度的差異以及中心像素強度組成，分別由f1和f2濾波器得出。
        self.__f1 = np.array([[1,1,1],
                              [1,-8,1],
                              [1,1,1]])
        self.__f2 = np.array([[0,0,0],
                              [0,1,0],
                              [0,0,0]])
        #方向反映局部窗內灰度變化的空間分布信息。通過局部窗內水平方向與垂直方向上鄰域像素點
        #的灰度差值比值的反正切變換來描述。 方向分為豎直方向和水平方向，由f3和f4濾波器得出。
        self.__f3 = np.array([[0,-1,0],
                              [0,0,0],
                              [0,1,0]])
        self.__f4 = np.array([[0,0,0],
                              [1,0,-1],
                              [0,0,0]])
        #方向量化個數   如果修改，也需要修改__classify_fai函數
        self.__T = 12
        #差分激勵量化個數   如果修改，也需要修改__classify_epc函數
        self.__M= 8 
        #每個頻段上將差分激勵均勻地量化為S格  可以修改
        self.__S = 4
    
    def calc_hist(self,image,concatenate=True):
        '''
        輸出統計直方圖 形狀[M,T,S]或[MxTxS,]
        
        args:
            image：輸入灰度圖
            concatenate：表明輸出直方圖是否合並為一維
        '''        
        his = np.zeros((self.__M,self.__T,self.__S),np.float32)
        rows,cols = image.shape[:2]
        #計算像素點個數
        sum_pix = rows*cols
        epc,theta = self.__compute(image)
        for i in range(rows):
            for j in range(cols):
                #把差分激勵分為8個區間，分別對應一個類別
                m = self.__classify_epc(epc[i][j])
                t = theta[i][j]
                s = self.__classify_s(epc[i][j],m)
                his[m-1][t-1][s-1] += 1 
        #歸一化
        his /= sum_pix
        if concatenate:
            his = np.reshape(his,-1)            
        return his

        
    def __compute(self,image):
        '''
        計算每個像素點的差分激勵和方向
        '''
        rows,cols = image.shape[:2]
        #用於保存每個像素點對應的差分激勵算子
        epc = np.zeros((rows,cols),dtype=np.float32)
        #用於保存每個像素點對應的方向算子
        theta = np.zeros((rows,cols),dtype=np.float32)
        
        '''
        計算差分激勵ξ
        '''
        v1 = cv2.filter2D(image,cv2.CV_16SC1,self.__f1)    
        v2 = cv2.filter2D(image,cv2.CV_16SC1,self.__f2)          
        for i in range(rows):
            for j in range(cols):
                #-π/2~π/2
                epc[i][j] = np.arctan(v1[i][j]/(v2[i][j]+0.0001))
                     
        '''
        計算每個像素點的方向Φ，把方向[0,2pi]均勻地量化為12個區間
        '''
        v3 = cv2.filter2D(image,cv2.CV_16SC1,self.__f3)
        v4 = cv2.filter2D(image,cv2.CV_16SC1,self.__f4)          
        for i in range(rows):
            for j in range(cols):
                theta[i][j] = np.arctan(v3[i][j]/((v4[i][j])+0.0001))
                if v3[i][j]>0 and v4[i][j]>0:
                    pass
                elif v3[i][j]>0 and v4[i][j]<0:
                    theta[i][j] = theta[i][j]+2*np.pi
                else:
                    theta[i][j]=theta[i][j]+np.pi
                theta[i][j] = self.__classify_fai(theta[i][j])             
        theta = theta.astype(np.uint8)
        return epc,theta
    
    def __classify_fai(self,value):
        '''
        把方向值value分類  類別為1、2、3、4、5、6、7、8、9、10、11、12 
        
        args:
            value：數值 0~2π之間
        '''
        if value >= 0 and value < 0.15*np.pi:            
            return 1
        elif value >= np.pi*0.15 and value < 0.35*np.pi:
            return 2
        elif value >= np.pi*0.35 and value < 0.5*np.pi:
            return 3
        elif value >= np.pi*0.5 and value < 0.65*np.pi:
            return 4
        elif value >= 0.65*np.pi and value < 0.85*np.pi:
            return 5
        elif value >= np.pi*0.85 and value < np.pi:
            return 6
        elif value >= np.pi and value < 1.15*np.pi:
            return 7
        elif value >= 1.15*np.pi and value < 1.35*np.pi:
            return 8
        elif value >= np.pi*1.35 and value < 1.5*np.pi:
            return 9
        elif value >= 1.5*np.pi and value < 1.65*np.pi:
            return 10
        elif value >= 1.65*np.pi and value < 1.85*np.pi:
            return 11
        else:
            return 12

    def __classify_epc(self,value):
        '''
        把差分激勵值value分類，划分為8個區間  類別為1、2、3、4、5、6、7、8 
        
        args:
            value：數值 -π/2~π/2之間
        '''
        if value >= np.pi*(-0.5) and value < (-0.3)*np.pi:
            return 1
        elif value >= np.pi*(-0.3) and value < (-0.15)*np.pi:
            return 2
        elif value >= np.pi*(-0.15) and value < (-0.05)*np.pi:
            return 3
        elif value >= np.pi*(-0.05) and value < 0:
            return 4
        elif value >= 0 and value < 0.05*np.pi:
            return 5
        elif value >= np.pi*0.05 and value < 0.15*np.pi:
            return 6
        elif value >= np.pi*0.15 and value < 0.3*np.pi:
            return 7
        else:
            return 8

    def __classify_s(self,value,label):
        '''
        將每個區間的差分激勵再次划分為S格
        
        args：
            value：差分激勵值   -π/2~π/2之間
            label：當前所屬區間  1,2,3,...,8
        '''
        if label == 1:
            space = ((-0.3)*np.pi - (-0.5)*np.pi)/self.__S
            return int((value - (-0.5)*np.pi)/space)+1
        elif label == 2:
            space = ((-0.15)*np.pi - (-0.3)*np.pi)/self.__S
            return int((value - (-0.3)*np.pi)/space)+1
        elif label == 3:
            space = ((-0.05)*np.pi - (-0.15)*np.pi)/self.__S
            return int((value - (-0.15)*np.pi)/space)+1
        elif label == 4:
            space = 0-(-0.05)*np.pi/self.__S
            return int((value - (-0.05)*np.pi)/space)+1
        elif label == 5:
            space = 0.05*np.pi /self.__S
            return int(value /space)+1
        elif label == 6:
            space = (0.15*np.pi - 0.05*np.pi)/self.__S
            return int((value - 0.05*np.pi)/space)+1
        elif label == 7:
            space = (0.3*np.pi - 0.15*np.pi)/self.__S
            return int((value - 0.15*np.pi)/space)+1
        else:
            space = (0.5*np.pi - 0.3*np.pi)/self.__S
            n = int((value - 0.3*np.pi)/space)+1
            if n == self.__S+1:
                n = self.__S
            return n
            
   
def  draw_hist(hist):
    '''
    繪制直方圖
    '''
    #首先先創建一個黑底的圖像，為了可以顯示彩色，所以該繪制圖像是一個8位的3通道圖像  
    #圖像高為100，寬為直方圖的長度
    width = len(hist)    
    height = 100
    draw_image = np.zeros((height,width,3),dtype= np.uint8)
    #獲取最大值
    max_value = np.max(hist)    
    #數值量化
    value = np.asarray((hist*0.9/max_value*height),dtype=np.int8)
    for i in range(width):
        cv2.line(draw_image,(i,height-1),(i,height-1-value[i]),(255,0,0))
    #顯示
    cv2.imshow('hist',draw_image)
    cv2.waitKey(0)
    cv2.destroyAllWindows()
    
    
if __name__=='__main__':
    image = cv2.imread('./image/match1.jpg')
    image = cv2.resize(image,dsize=(600,400))
    imgray = cv2.cvtColor(image, cv2.COLOR_BGR2GRAY)
    wld = WLD()
    hist = wld.calc_hist(imgray)
    print(hist.shape)
    print(hist)
    print(np.sum(hist))
    draw_hist(hist)