圖像特征及提取

• • 1. 傳統圖像特征
  • • 1.1 統計特征---穿透數、不變矩陣、粗網格、方向線素
    • 1.2 紋理特征---LBP, Gabor, HOG, GLCM
  • 2. 圖像深度特征

1. 傳統圖像特征

1.1 統計特征—穿透數、不變矩陣、粗網格、方向線素

筆畫穿透數（stroke penetration）是一種全局統計特征（global statistical feature），它描述字元的各部分筆劃的疏密程度，能夠提供比較完整的字元資訊。在網格字元圖像中，從不同方向進行掃描。根據輪廓交叉(silhouette intersection)的數量，定義穿透數特征函數H。選取水準方向和垂直方向，在每個方向上選取n(n=16)個特征值，然後就可擷取到 2 n 2n 2n個特征向量函數 H ， H H，H H，H由以下公式定義：

H = ( h 11 , h 12 , . . . , h 1 n , h 21 , h 22 , . . . , h 2 n ) , n = ( 1 , 2 , . . . , 16 ) − − − − ( 1 ) H=(h_{11},h_{12},...,h_{1n},h_{21},h_{22},...,h_{2n}), n=(1,2,...,16) ----(1) H=(h11​,h12​,...,h1n​,h21​,h22​,...,h2n​),n=(1,2,...,16)−−−−(1)

假設字元的水準和垂直坐标被規範化為 1 ∼ p 1∼p 1∼p，設H函數的每個初值為 p n = p / n p_n=p/n pn​=p/n ，則集合H中的 h n h_n hn​ 表示為：

h n = n / p ∗ ∑ t = 1 p n h ( k + p s − 1 ) − − − − ( 2 ) h_n=n/p*\sum_{t=1}^{p_n}h(k+p_s-1)----(2) hn​=n/p∗t=1∑pn​​h(k+ps​−1)−−−−(2)

其中 s = 1 , 2 , t = 1 , 2 , . . . , n s=1,2, t=1,2,...,n s=1,2,t=1,2,...,n， h n h_n hn​為字元穿透數特征向量。

不變矩陣（moment）是全局統計特征，其在圖像檢索（image retrieval）和圖像識别（image recognition）中得到了廣泛的應用。由于圖像區域的某些矩值對于平移（translation）、旋轉（rotate）、縮放(scale)等幾何操作是不可改變的，是以矩值的表示方法在模式識别中非常重要。對于二維連續函數（2D continuous function） f ( x , y ) ， f(x,y)， f(x,y)，階矩（order moment）(j+k)定義如下：

m j k = ∫ − ∞ + ∞ ∫ − ∞ + ∞ x j y k f ( x , y ) d x d y − − − − ( 3 ) m_{jk}=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} x^jy^kf(x,y)d_xd_y ----(3) mjk​=∫−∞+∞​∫−∞+∞​xjykf(x,y)dx​dy​−−−−(3)

其中 j 和 k j和k j和k都是非負整數，是以得到了一個無限的矩集。為了描述字元圖像的形狀，在函數f(x,y)中假設字元的函數值為1，背景函數值為0。是以該函數隻反映字元對象的形狀，忽略灰階級細節。參數j+k稱為矩序。零階矩（zeroth order moment）則是字元對象的面積，表示為

m 00 = ∫ − ∞ + ∞ ∫ − ∞ + ∞ f ( x , y ) d x d y − − − − ( 4 ) m_{00}=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} f(x,y)d_xd_y ----(4) m00​=∫−∞+∞​∫−∞+∞​f(x,y)dx​dy​−−−−(4)

當滿足條件 j = 1 ， k = 0 時 j=1，k=0時 j=1，k=0時， m 10 m_{10} m10​為二值字元圖像中所有點的水準坐标x之和。同樣， m 01 m_{01} m01​ 是圖像中所有點的垂直坐标y之和。

設 ( x ‾ , y ‾ ) ( x ‾ = m 10 m 00 , y ‾ = m 01 m 00 ) (\overline x,\overline y)(\overline x=\frac {m_{10}} {m_{00}},\overline y= \frac {m_{01}} {m_{00}}) (x,y​)(x=m00​m10​​,y​=m00​m01​​)為二值字元圖像的質心坐标，則中心矩（central moment）定義為

μ j k = ∫ − ∞ + ∞ ∫ − ∞ + ∞ ( x − x ‾ ) j ( y − y ‾ ) k f ( x , y ) d x d y − − − − ( 5 ) \mu_{jk}=\int_{-\infty}^{+\infty} \int_{-\infty}^{+\infty} {(x-\overline x)^j}{(y-\overline y)^k}f(x,y)d_xd_y ----(5) μjk​=∫−∞+∞​∫−∞+∞​(x−x)j(y−y​)kf(x,y)dx​dy​−−−−(5)

粗網格（grid pixels）特征是局部統計特征，其反映了圖像形狀的分布情況。粗網格統計值即為字元圖像中每個網格中黑色像素所占的百分比，每個網格分别反映字元的某一部分特征。對樣本進行預處理後，将字元劃分為8*8的網格圖像，每個坐标的像素點 g ( i , j ) g(i,j) g(i,j)可以表示為 g ( i , j ) = { 1 , b l a c k 0 , w h i t e } − − − − ( 6 ) g(i,j)=\{_{1,black} ^{0, white}\}----(6) g(i,j)={1,black0,white​}−−−−(6)

對每個網格内的前景像素進行統計，可得到64維粗網格特征向量。

方向線素（direction line element）是一種局部統計特征（local statistical feature），即輪廓分解後的網格方向特征，在中文和英文等識别中表現的效果較好。方向線素特征描述筆畫在不同方向上的貢獻度，筆畫方向量化為如下圖左側所示的4個方向，分别是垂直，水準和傾斜正負 4 5 o 45^o 45o。

對于3*3網格中的中心黑色像素，可以觀察到兩種情況：一種是配置設定了一種線元素(下圖(a)-圖(d))；另一種是三個黑色像素相連，但配置設定了兩種類型的線元素(圖(e)-圖(l))。例如，在圖(f)中，同時配置設定 4 5 o 45^o 45o的線元素和垂直線元素，以此量化不同方向的線元素。

将字元圖像劃分為具有16×16像素的7×7子區域，然後将每個子區域分成4個嵌套塊，表示為A、B、C和D，尺寸分别為16×16、12×12、8×8和4×4。對于每個子區域，定義一個四維向量 X = ( x 1 , x 2 , x 3 , x 4 ) T X=(x_1,x_2,x_3,x_4)^T X=(x1​,x2​,x3​,x4​)T ，其中 x 1 x_1 x1​、 x 2 x_2 x2​ 、 x 3 x_3 x3​ 和 x 4 x_4 x4​分别表示四個方向上的方向線素的數量，四個方向包括 0 o 、 4 5 o 、 9 0 o 和 13 5 o 。 0^o、45^o、90^o和135^o。 0o、45o、90o和135o。

為了避免受到字元圖像位置改變的影響，将對每個區域配置設定權重。對于16×16像素的子區域，每個子區域配置設定的權重可表示為

其中， x j A x_j ^A xjA​ 、 x j B x_j ^B xjB​ 、 x j C x_j ^C xjC​ 、 x j D x_j ^D xjD​ 分别表示區域A、區域B、區域C、區域D中某個方向上的方向線素的數量。每個子區域可以得到一個四維特征向量，每個字元圖像得到7×7×4=196維方向線素特征向量。

1.2 紋理特征—LBP, Gabor, HOG, GLCM

Local Binary Pattern(LBP), Gabor, Histogram of Oriented Gradient(HOG), 灰階共生矩陣(Gray Level CO-Occurrence Matrix, GLCM)

LBP：可以看之前的這篇部落格，以及本篇的代碼公開在GitHub。

Gabor : 1946年提出Gabor濾波器，可以在頻域的不同尺度、不同方向上提取相關的特征。對于圖像的邊緣敏感，能夠提供良好的方向選擇和尺度選擇特性，Gabor濾波器的這一性質，使得其在視覺領域中經常被用來作圖像的預處理。

Gabor濾波器示意圖，4種角度6種方向：

以及生成的24張濾波圖像。

#建構Gabor濾波器
def build_filters():
     filters = []
     ksize = [9,10,11,12,13] # gabor尺度，6個
     lamda = np.pi/6.0         # 波長
     for theta in np.arange(0, np.pi, np.pi / 8): #gabor方向，0°，45°，90°，135°，共四個
         for K in range(5):
             kern = cv2.getGaborKernel((ksize[K], ksize[K]), 7.0, theta, lamda, 0.5, 0, ktype=cv2.CV_32F)
             kern /= 1.5*kern.sum()
             filters.append(kern)
     return filters
           

HOG：2005年 CVPR 行人識别。思想為局部目标的表象和形狀能夠被梯度或邊緣的方向密度分布很好地描述，統計邊緣梯度資訊。

步驟：1.預處理 gama矯正 灰階化等

2. 計算梯度 計算圖像橫縱(取絕對值)方向的梯度，并據此得出每個像素位置的梯度方向 θ \theta θ及梯度強度 g g g。

θ = a r c t a n g x g y g = g x 2 + g y 2 \theta=arctan \frac {g_x} {g_y} \qquad g=\sqrt{g_x ^2 +g_y ^2} θ=arctangy​gx​​g=gx2​+gy2​

​

3. 計算梯度直方圖 按照梯度的方向和強度，将cell中對應的像素劃分至0~180度的9個bins中，得到9維的向量。

4. 向上合并 四個cell進行歸一化

HOG Python實作：

from skimage.feature import hog, local_binary_pattern
import cv2

img = cv2.imread('/1.jpg', cv2.IMREAD_GRAYSCALE)
img = cv2.resize(img, (512, 512))
# gray = rgb2gray(image) / 255.0
fd = hog(img, orientations=12, block_norm='L1', pixels_per_cell=[10, 10], cells_per_block=[4, 4], visualize=False, transform_sqrt=True)
print(fd.shape)
radius = 3
n_point = radius * 16
lbp = local_binary_pattern(img,n_point,radius,'default')
cv2.imshow('img', lbp)
cv2.waitKey(0)
           

import cv2
import numpy as np
import math
import matplotlib.pyplot as plt

class Hog_descriptor():
    def __init__(self, img, cell_size=16, bin_size=8):
        self.img = img
        self.img = np.sqrt(img / np.max(img))
        self.img = img * 255
        self.cell_size = cell_size
        self.bin_size = bin_size
        self.angle_unit = 360 // self.bin_size
        assert type(self.bin_size) == int, "bin_size should be integer,"
        assert type(self.cell_size) == int, "cell_size should be integer,"
        assert type(self.angle_unit) == int, "angle_unit should be divisible by 360"

    def extract(self):
        height, width = self.img.shape
        print(height / self.cell_size)
        gradient_magnitude, gradient_angle = self.global_gradient()
        gradient_magnitude = abs(gradient_magnitude)
        cell_gradient_vector = np.zeros((height // self.cell_size, width // self.cell_size, self.bin_size))
        for i in range(cell_gradient_vector.shape[0]):
            for j in range(cell_gradient_vector.shape[1]):
                cell_magnitude = gradient_magnitude[i * self.cell_size:(i + 1) * self.cell_size,
                                 j * self.cell_size:(j + 1) * self.cell_size]
                cell_angle = gradient_angle[i * self.cell_size:(i + 1) * self.cell_size,
                             j * self.cell_size:(j + 1) * self.cell_size]
                cell_gradient_vector[i][j] = self.cell_gradient(cell_magnitude, cell_angle)

        hog_image = self.render_gradient(np.zeros([height, width]), cell_gradient_vector)
        hog_vector = []
        for i in range(cell_gradient_vector.shape[0] - 1):
            for j in range(cell_gradient_vector.shape[1] - 1):
                block_vector = []
                block_vector.extend(cell_gradient_vector[i][j])
                block_vector.extend(cell_gradient_vector[i][j + 1])
                block_vector.extend(cell_gradient_vector[i + 1][j])
                block_vector.extend(cell_gradient_vector[i + 1][j + 1])
                mag = lambda vector: math.sqrt(sum(i ** 2 for i in vector))
                magnitude = mag(block_vector)
                if magnitude != 0:
                    normalize = lambda block_vector, magnitude: [element / magnitude for element in block_vector]
                    block_vector = normalize(block_vector, magnitude)
                hog_vector.append(block_vector)
        return hog_vector, hog_image

    def global_gradient(self):
        gradient_values_x = cv2.Sobel(self.img, cv2.CV_64F, 1, 0, ksize=5)
        gradient_values_y = cv2.Sobel(self.img, cv2.CV_64F, 0, 1, ksize=5)
        gradient_magnitude = cv2.addWeighted(gradient_values_x, 0.5, gradient_values_y, 0.5, 0)
        gradient_angle = cv2.phase(gradient_values_x, gradient_values_y, angleInDegrees=True)
        return gradient_magnitude, gradient_angle

    def cell_gradient(self, cell_magnitude, cell_angle):
        orientation_centers = [0] * self.bin_size
        for i in range(cell_magnitude.shape[0]):
            for j in range(cell_magnitude.shape[1]):
                gradient_strength = cell_magnitude[i][j]
                gradient_angle = cell_angle[i][j]
                min_angle, max_angle, mod = self.get_closest_bins(gradient_angle)
                orientation_centers[min_angle] += (gradient_strength * (1 - (mod / self.angle_unit)))
                orientation_centers[max_angle] += (gradient_strength * (mod / self.angle_unit))
        return orientation_centers

    def get_closest_bins(self, gradient_angle):
        idx = int(gradient_angle / self.angle_unit)
        mod = gradient_angle % self.angle_unit
        return idx, (idx + 1) % self.bin_size, mod

    def render_gradient(self, image, cell_gradient):
        cell_width = self.cell_size / 2
        max_mag = np.array(cell_gradient).max()
        for x in range(cell_gradient.shape[0]):
            for y in range(cell_gradient.shape[1]):
                cell_grad = cell_gradient[x][y]
                cell_grad /= max_mag
                angle = 0
                angle_gap = self.angle_unit
                for magnitude in cell_grad:
                    angle_radian = math.radians(angle)
                    x1 = int(x * self.cell_size + magnitude * cell_width * math.cos(angle_radian))
                    y1 = int(y * self.cell_size + magnitude * cell_width * math.sin(angle_radian))
                    x2 = int(x * self.cell_size - magnitude * cell_width * math.cos(angle_radian))
                    y2 = int(y * self.cell_size - magnitude * cell_width * math.sin(angle_radian))
                    cv2.line(image, (y1, x1), (y2, x2), int(255 * math.sqrt(magnitude)))
                    angle += angle_gap
        return image

img = cv2.imread('/1.jpg', cv2.IMREAD_GRAYSCALE)
img = cv2.resize(img, (512, 512))

hog = Hog_descriptor(img, cell_size=8, bin_size=8)
vector, image = hog.extract()
print(np.array(vector).shape)
plt.imshow(image, cmap=plt.cm.gray)
plt.show()
cv2.imshow('img', image)
cv2.waitKey(0)
           

GLCM：1973年Haralick等人提出，是一種通過研究灰階的空間相關特性來描述紋理的常用方法，由于紋理是由灰階分布在空間位置上反複出現而形成的，因而在圖像空間中相隔某距離的兩像素之間會存在一定的灰階關系，即圖像中灰階的空間相關特性。

定義：從灰階為 i 的像素點出發，距離(dx，dy)的另一像素點灰階為 j 的的機率。

灰階直方圖是對圖像上單個像素具有某個灰階進行統計的結果，而灰階共生矩陣是對圖像上保持某距離的兩像素分别具有某灰階的狀況進行統計得到的。

取圖像(N×N)中任意一點（x，y）及偏離它的另一點（x+a，y+b），設該點對的灰階值為 （g1，g2）。令點（x，y） 在圖像上移動，則會得到多個（g1，g2）值，設灰階值的級數為 k，則（g1，g2） 的組合共有 k ∗ k k*k k∗k種。對于整幅圖像，統計出每一種 （g1，g2）值出現的次數，然後排列成一個方陣，再用（g1，g2） 出現的總次數将它們歸一化為出現的機率 P ( g 1 ， g 2 ） P_{(g1，g2）} P(g1，g2）​ ，這樣的方陣稱為灰階共生矩陣。

其中，d：用像素數量表示的相對距離； θ \theta θ : 一般考慮四個方向，分别為 0°,45°,90°,135° ; i , j = 0 , 1 , 2 , . . . , L − 1 i, j =0,1,2,...,L-1 i,j=0,1,2,...,L−1 ;（x，y）為圖像中的像素坐标，L為圖像灰階級的數目。可以簡單了解為灰階圖像中某種形狀的像素對，在圖像中出現的次數。

1. 角二階矩(Angular Second Moment, ASM)，又稱能量，是圖像灰階分布均勻程度和紋理粗細的一個度量。

   

2. 熵(Entropy, ENT)，表示圖像中紋理的非均勻程度或複雜程度。熵是圖像所具有的資訊量的度量，紋理資訊也屬于圖像的資訊，是一個随機性的度量。當共生矩陣中所有元素有最大的随機性、空間共生矩陣中所有值幾乎相等時，共生矩陣中元素分散分布時，熵較大。

   

3. 反差分矩陣(Inverse Differential Moment, IDM)，反映圖像紋理的同質性，度量圖像紋理局部變化的多少。

   

   Python 代碼實作：

def get_inputs(input):  # s為圖像路徑
    # 得到共生矩陣，參數：圖像矩陣，距離，方向，灰階級别，是否對稱，是否标準化
    # glcm = greycomatrix(input, [16, 32, 48, 64, 80, 96, 112, 128, 144, 160], [0, np.pi / 4, np.pi / 2, np.pi * 3 / 4], 256, symmetric=True, normed=True)
    glcm = greycomatrix(input, [1], [0, np.pi / 4, np.pi / 2, np.pi * 3 / 4], 256, symmetric=True, normed=True)
    # 得到共生矩陣統計值，http://tonysyu.github.io/scikit-image/api/skimage.feature.html#skimage.feature.greycoprops
    #temp = greycoprops(glcm,'contrast')
    #temp = greycoprops(glcm,'dissimilarity')
    #temp = greycoprops(glcm, 'homogeneity')
    #temp = greycoprops(glcm, 'correlation')
    temp = greycoprops(glcm, 'energy')
    print(temp)
    #temp = greycoprops(glcm, 'ASM')
    # temp1 = np.array(temp)
    # temp2 = list(chain.from_iterable(temp1))
    return temp
           

if __name__ == '__main__':
    filters = build_filters()
    print(len(filters))
    a = 0 
    contrast = np.zeros((40000,40), dtype=np.float32)
    dissimilarity = np.zeros((40000,40), dtype=np.float32)
    homogeneity = np.zeros((40000,40), dtype=np.float32)
    ASM = np.zeros((40000,40), dtype=np.float32)

    for filename in os.listdir('D:/'):
        img = cv2.imread('D:/' + filename, 0)
        print(filename)
        for i in range(len(filters)):
            accum = np.zeros_like(img)
            for kern in filters[i]:
                fimg = cv2.filter2D(img, cv2.CV_8UC1, kern)
                accum = np.maximum(accum, fimg, accum)
                # print(i)
                # cv2.imwrite(str(i) + ".jpg", accum)
            contrast1, dissimilarity1, homogeneity1, ASM1 = get_inputs(accum)
            # print(contrast, np.mean(contrast))
            # energy.append(np.mean(temp))
            # for i in range(40):
            ASM[a, i] = np.mean(ASM1)
            contrast[a, i] = np.mean(contrast1)
            dissimilarity[a, i] = np.mean(dissimilarity1)
            homogeneity[a, i] = np.mean(homogeneity1)
        a += 1    
    contrast = np.array(contrast)
    dissimilarity = np.array(dissimilarity)
    homogeneity = np.array(homogeneity)
    ASM = np.array(ASM)
    # energy = list(chain.from_iterable(energy))
    # print('mena:', mean)
    contrast = pd.DataFrame(data=contrast)
    contrast.to_csv("./feature-/contrast.csv", encoding='utf-8', index=False)
    dissimilarity = pd.DataFrame(data=dissimilarity)
    dissimilarity.to_csv("./feature-/dissimilarity.csv", encoding='utf-8', index=False)
    homogeneity = pd.DataFrame(data=homogeneity)
    homogeneity.to_csv("./feature-/homogeneity.csv", encoding='utf-8', index=False)
    ASM = pd.DataFrame(data=ASM)
    ASM.to_csv("./feature-/ASM.csv", encoding='utf-8', index=False)
           

2. 圖像深度特征

設計神經網絡模型來挖掘圖像更深、更為抽象的特征。CNN随着深度的增加時所表達的進階語義資訊比較強，而淺層的卷積結構對物體的輪廓刻畫比較明顯

下圖，原圖經卷積後随深度增加，特征随之抽象。代碼公開在GitHub

model.save('./m_mnist.h5')
from keras.models import load_model
model = load_model('./m_mnist.h5')

def conv_output(model, layer_name, img):
    """Get the output of conv layer.

    Args:
           model: keras model.
           layer_name: name of layer in the model.
           img: processed input image.

    Returns:
           intermediate_output: feature map.
    """
    # this is the placeholder for the input images
    input_img = model.input
    print('shape', input_img.shape)
    try:
        # this is the placeholder for the conv output
        out_conv = model.get_layer(name=layer_name).output
    except:
        raise Exception('Not layer named {}!'.format(layer_name))

    # get the intermediate layer model
    intermediate_layer_model = Model(inputs=input_img, outputs=out_conv)
    # get the output of intermediate layer model
    intermediate_output = intermediate_layer_model.predict(img)

    return intermediate_output[0]

def vis_conv(images, n, name, t):
    """visualize conv output and conv filter.
    Args:
           img: original image.
           n: number of col and row.
           t: vis type.
           name: save name.
    """
    size = 64
    margin = 5

    if t == 'filter':
        results = np.zeros((n * size + 7 * margin, n * size + 7 * margin, 3))
    if t == 'conv':
        results = np.zeros((n * size + 7 * margin, n * size + 7 * margin))

    for i in range(n):
        for j in range(n):
            if t == 'filter':
                filter_img = images[i + (j * n)]
            if t == 'conv':
                filter_img = images[..., i + (j * n)]
            filter_img = cv2.resize(filter_img, (size, size))

            # Put the result in the square `(i, j)` of the results grid
            horizontal_start = i * size + i * margin
            horizontal_end = horizontal_start + size
            vertical_start = j * size + j * margin
            vertical_end = vertical_start + size
            if t == 'filter':
                results[horizontal_start: horizontal_end, vertical_start: vertical_end, :] = filter_img
            if t == 'conv':
                results[horizontal_start: horizontal_end, vertical_start: vertical_end] = filter_img

    # Display the results grid
    plt.imshow(results)
    # plt.savefig('images\{}_{}.jpg'.format(t, name), dpi=600)
    plt.show()
    
def visualize_feature_map(img_batch):
    feature_map = img_batch
    print(feature_map.shape)
    feature_map_combination = []
    plt.figure()
 
    num_pic = feature_map.shape[2]
    row, col = get_row_col(num_pic)
 
    for i in range(0, num_pic):
        feature_map_split = feature_map[:, :, i]
        feature_map_combination.append(feature_map_split)
        plt.subplot(row, col, i + 1)
        plt.imshow(feature_map_split)
        axis('off')
    plt.savefig('./mu_feature_map.png')
 
    # 各個特征圖按1：1 疊加
    feature_map_sum = sum(ele for ele in feature_map_combination)
    plt.imshow(feature_map_sum)
    plt.savefig("./feature_map_sum.jpg")
    plt.show()
# model 注意替換
img = X_test[0].reshape(1, img_rows, img_cols, 1)
img = conv_output(model, 'conv2d_8', img)
vis_conv(img, 8, 'conv2d_8', 'conv')

from keras.models import load_model
from keras.utils import CustomObjectScope
with CustomObjectScope({'SpatialPyramidPooling': SpatialPyramidPooling}):
    model = load_model('./m_repvgg40e.h5')
model = Model(inputs=model.input, outputs=model.get_layer('add2').output)    
img1 = cv2.imread('D:/guji/yi69_00200_22.jpg', 0)
block_pool_features = model.predict(img1)
print(block_pool_features.shape)
feature = block_pool_features.reshape(block_pool_features.shape[1:])
visualize_feature_map(feature)