画素毎にテクスチャの特徴抽出をするのに、全ての画素に対してその周辺窓のGLCMを計算したい。この計算をするのに、pythonなのに画素の二重forループ使って計算する方法があるけど、遅すぎて使い物にならない。そこで、numpyを使って画素毎のGLCM計算をスクラッチから書いてみて、scikit-imageのものと比較してみる。

まずは結果から

f:id:tzmi:20200206215059j:plain

scikit-imageではこのうち、mean, std, energyの計算ができない。

次に計算速度の違い。numpyにしたら300倍早くなった。

153.572 sec  #for loopの場合
0.505 sec # numpyの場合

GLCMって何？

GLCMはGray Level Co-occurrence Matrices の略で元論文は R. M. Hararick et al. 1973, "Textural Features for Image Classification" と非常に古い。ある画像領域のテクスチャを示す行列(GLCM)を抽出し、この行列から様々な特徴量を計算する。GLCMはある画像領域での隣接画素との関係から得られる２次元ヒストグラムである。

計算を高速化するためには各画素のビット数を落として隣接画素との大小関係をおおざっぱに比べる。例えば0~255の8bitを0-７の3bitくらいまで落として計算する。テクスチャだけでなく明るい領域や暗い領域を区別できる特徴も計算できる。

技術の詳細な解説記事は例えば以下。 qiita.com

この記事では画像全体ではなく、画素毎のGLCM計算を高速化する方法について紹介する。

numpyを使った画素毎GLCM計算の高速化

ということで、numpyバージョンを実装していく。方針としてはim2colとかと同じでforループしないといけない部分を最小限に抑えながら計算する。画素毎ではなくmatrixの成分毎にループをまわすだけ（比較画素のビット数を落とすことが前提）。scikit-imageには方向というパラメータがあるが、上か右かくらいしか使わないし方向変えても写真のような画像ではあんまり変わらないから右の画素と比較する。距離も１画素とするが、距離をパラメータ化したかったら、opencvとかで画像サイズを変えれば実現できると思う。

注意点として、numpyを使うときはメモリ容量が大きいPCで実行することをお勧めする。メモリ容量が小さいとnumpyでGLCMの配列を確保するところで詰む。下記のコードならlevel=８なので16MBくらいで大丈夫だけど、level=128とかするとメモリが4GB必要になる。levelの値は大きくしないように注意。

import numpy as np
import matplotlib.pyplot as plt
from skimage import data
import cv2

import matplotlib as mpl
from time import time

def main():
    mpl.rc('image', cmap='inferno')
    pass


if __name__ == '__main__':
    main()

    kernel_size = 5
    levels = 8
    symmetric = False
    normed = True

    img = data.camera()

    # binarize
    img_bin = img//(256//levels) # [0:255]->[0:7]

    # calc_glcm                                                                 
    st = time()
    h,w = img.shape

    glcm = np.zeros((h,w,levels,levels), dtype=np.uint8)
    kernel = np.ones((kernel_size, kernel_size), np.uint8)
    img_bin_r = np.append(img_bin[:,1:], img_bin[:,-1:], axis=1)
    for i in range(levels):
        for j in range(levels):
            mask = (img_bin==i) & (img_bin_r==j)
            mask = mask.astype(np.uint8)
            glcm[:,:,i,j] = cv2.filter2D(mask, -1, kernel)

    glcm = glcm.astype(np.float32)

    if symmetric:
        glcm += glcm[:,:,::-1, ::-1]

    if normed:
        glcm = glcm/glcm.sum(axis=(2,3), keepdims=True)

    # martrix axis
    axis = np.arange(levels, dtype=np.float32)+1
    w = axis.reshape(1,1,-1,1)
    x = np.repeat(axis.reshape(1,-1), levels, axis=0)
    y = np.repeat(axis.reshape(-1,1), levels, axis=1)

    # GLCM contrast                                                             
    glcm_contrast = np.sum(glcm*(x-y)**2, axis=(2,3))

    # GLCM dissimilarity                                                        
    glcm_dissimilarity = np.sum(glcm*np.abs(x-y), axis=(2,3))

    # GLCM homogeneity                                                          
    glcm_homogeneity = np.sum(glcm/(1.0+(x-y)**2), axis=(2,3))

    # GLCM energy & ASM                                                         
    glcm_asm = np.sum(glcm**2, axis=(2,3))

    # GLCM entropy                                                              
    ks = 5 # kernel_size                                                        
    pnorm = glcm / np.sum(glcm, axis=(2,3), keepdims=True) + 1./ks**2
    glcm_entropy = np.sum(-pnorm * np.log(pnorm), axis=(2,3))

    print('%.3f sec' %(time()-st))


    # GLCM mean                                                                 
    glcm_mean = np.mean(glcm*w, axis=(2,3))

    # GLCM std                                                                  
    glcm_std = np.std(glcm*w, axis=(2,3))

    # GLCM energy                                                               
    glcm_energy = np.sqrt(glcm_asm)

    # GLCM max                                                                  
    glcm_max = np.max(glcm, axis=(2,3))


    # plot                                                                      
    plt.figure(figsize=(10,4.5))
    fs = 15

    plotid = [1,2,3,4,5,6,7,8,9,10]
    outs =[img, glcm_mean, glcm_std,
           glcm_contrast, glcm_dissimilarity, glcm_homogeneity,
           glcm_asm, glcm_energy, glcm_max,
           glcm_entropy]
    titles = ['original', 'mean', 'std',
              'contrast', 'dissimilarity', 'homogeneity',
              'ASM', 'energy', 'max',
              'entropy']

    for i in range(10):
        plt.subplot(2,5,plotid[i])
        plt.tick_params(labelbottom=False, labelleft=False)
        plt.xticks([])
        plt.yticks([])
        plt.imshow(outs[i])
        plt.title(titles[i], fontsize=fs)

    plt.tight_layout(pad=0.5)
    plt.savefig('numpy_output.png')
    plt.show()