GLCM計算のためのPythonのスライディングウィンドウ

Question

GLCMアルゴリズムを使用して衛星画像でテクスチャ解析をしようとしています。シキット画像のドキュメントは非常に役立ちますが、GLCMの計算には、画像上にウィンドウサイズをループする必要があります。これはPythonでは遅すぎます。私はスライディングウインドウについてのstackoverflowに関する多くの記事を見つけましたが、計算はこれまでどおりです。私は下に示す例を持っている、それは動作しますが、永遠にかかる。私はこれは私が私のために良い解決策かもしれませんskimage.util.view_as_windows方法で出くわした

image = np.pad(image, int(win/2), mode='reflect') 
row, cols = image.shape 
feature_map = np.zeros((M, N)) 

for m in xrange(0, row): 
    for n in xrange(0, cols): 
     window = image[m:m+win, n:n+win] 
     glcm = greycomatrix(window, d, theta, levels) 
     contrast = greycoprops(glcm, 'contrast') 
     feature_map[m,n] = contrast 

をそれを行うの素朴な方法でなければならないと思います。私の問題は、私はGLCMを計算しようとすると、私が言うエラーが出る、ということである：

ValueError: The parameter image must be a 2-dimensional array

GLCMの画像の結果は、4D寸法とscikit-画像view_as_windows方法が唯一の2D配列を受け入れるを持っているためです。ここに私の試みです

win_w=40 
win_h=40 

features = np.zeros(image.shape, dtype='uint8') 
target = features[win_h//2:-win_h//2+1, win_w//2:-win_w//2+1] 
windowed = view_as_windows(image, (win_h, win_w)) 

GLCM = greycomatrix(windowed, [1], [0, np.pi/4, np.pi/2, 3*np.pi/4], symmetric=True, normed=True) 
haralick = greycoprops(GLCM, 'ASM') 

誰もが私はskimage.util.view_as_windowsメソッドを使用してGLCMを計算する方法についてのアイデアを持っていますか？

Answer 1

実行しようとしている特徴抽出は、コンピュータを大量に使用する作業です。スライディングウインドウの重なり合った位置で何回も共起マップを繰り返し計算するのではなく、共起マップを画像全体に対して1回だけ計算することによって、あなたの方法をスピードアップしました。

共起マップは、元の画像と同じサイズの画像のスタックであり、各ピクセルについて、強度レベルは、2つの強度の同時発生を符号化する整数で置き換えられる。すなわち、Ii atそのピクセルとオフセットピクセルでIj。共起マップは、オフセットを考慮した数の層（すなわち、すべての可能な距離 - 角度の組）を有する。共起マップを保持することにより、以前に計算された共起マップを再利用して各距離の隣接行列（GLCM）を得ることができるため、スライディングウィンドウの各位置でGLCMをゼロから計算する必要はありません角のペア。この方法は、スピードを大幅に向上させます。私が思いついた

溶液は、以下の機能に依存して、次の結果がランドサット画像から(250, 200)画素作物に対応

import numpy as np 
from skimage import io 
from scipy import stats 
from skimage.feature import greycoprops 

def offset(length, angle): 
    """Return the offset in pixels for a given length and angle""" 
    dv = length * np.sign(-np.sin(angle)).astype(np.int32) 
    dh = length * np.sign(np.cos(angle)).astype(np.int32) 
    return dv, dh 

def crop(img, center, win): 
    """Return a square crop of img centered at center (side = 2*win + 1)""" 
    row, col = center 
    side = 2*win + 1 
    first_row = row - win 
    first_col = col - win 
    last_row = first_row + side  
    last_col = first_col + side 
    return img[first_row: last_row, first_col: last_col] 

def cooc_maps(img, center, win, d=[1], theta=[0], levels=256): 
    """ 
    Return a set of co-occurrence maps for different d and theta in a square 
    crop centered at center (side = 2*w + 1) 
    """ 
    shape = (2*win + 1, 2*win + 1, len(d), len(theta)) 
    cooc = np.zeros(shape=shape, dtype=np.int32) 
    row, col = center 
    Ii = crop(img, (row, col), win) 
    for d_index, length in enumerate(d): 
     for a_index, angle in enumerate(theta): 
      dv, dh = offset(length, angle) 
      Ij = crop(img, center=(row + dv, col + dh), win=win) 
      cooc[:, :, d_index, a_index] = encode_cooccurrence(Ii, Ij, levels) 
    return cooc 

def encode_cooccurrence(x, y, levels=256): 
    """Return the code corresponding to co-occurrence of intensities x and y""" 
    return x*levels + y 

def decode_cooccurrence(code, levels=256): 
    """Return the intensities x, y corresponding to code""" 
    return code//levels, np.mod(code, levels)  

def compute_glcms(cooccurrence_maps, levels=256): 
    """Compute the cooccurrence frequencies of the cooccurrence maps""" 
    Nr, Na = cooccurrence_maps.shape[2:] 
    glcms = np.zeros(shape=(levels, levels, Nr, Na), dtype=np.float64) 
    for r in range(Nr): 
     for a in range(Na): 
      table = stats.itemfreq(cooccurrence_maps[:, :, r, a]) 
      codes = table[:, 0] 
      freqs = table[:, 1]/float(table[:, 1].sum()) 
      i, j = decode_cooccurrence(codes, levels=levels) 
      glcms[i, j, r, a] = freqs 
    return glcms 

def compute_props(glcms, props=('contrast',)): 
    """Return a feature vector corresponding to a set of GLCM""" 
    Nr, Na = glcms.shape[2:] 
    features = np.zeros(shape=(Nr, Na, len(props))) 
    for index, prop_name in enumerate(props): 
     features[:, :, index] = greycoprops(glcms, prop_name) 
    return features.ravel() 

def haralick_features(img, win, d, theta, levels, props): 
    """Return a map of Haralick features (one feature vector per pixel)""" 
    rows, cols = img.shape 
    margin = win + max(d) 
    arr = np.pad(img, margin, mode='reflect') 
    n_features = len(d) * len(theta) * len(props) 
    feature_map = np.zeros(shape=(rows, cols, n_features), dtype=np.float64) 
    for m in xrange(rows): 
     for n in xrange(cols): 
      coocs = cooc_maps(arr, (m + margin, n + margin), win, d, theta, levels) 
      glcms = compute_glcms(coocs, levels) 
      feature_map[m, n, :] = compute_props(glcms, props) 
    return feature_map 

---

DEMO

。私は2つの距離、4つの角度、2つのGLCM特性を考慮しました。これにより、各ピクセルについて16次元の特徴ベクトルが得られる。スライディングウィンドウが2乗され、その側面が2*win + 1ピクセル（このテストでは、win = 19の値が使用されたことに注意してください）。このサンプルの実行には、我々は、おそらく高い次元の配列をサポートするためにview_as_windowsを展開する必要があります;-)「永遠に」よりもかなり短くなっている、

In [331]: img.shape 
Out[331]: (250L, 200L) 

In [332]: img.dtype 
Out[332]: dtype('uint8') 

In [333]: d = (1, 2) 

In [334]: theta = (0, np.pi/4, np.pi/2, 3*np.pi/4) 

In [335]: props = ('contrast', 'homogeneity') 

In [336]: levels = 256 

In [337]: win = 19 

In [338]: %time feature_map = haralick_features(img, win, d, theta, levels, props) 
Wall time: 5min 53s  

In [339]: feature_map.shape 
Out[339]: (250L, 200L, 16L) 

In [340]: feature_map[0, 0, :]  
Out[340]: 
array([ 10.3314, 0.3477, 25.1499, 0.2738, 25.1499, 0.2738, 
     25.1499, 0.2738, 23.5043, 0.2755, 43.5523, 0.1882, 
     43.5523, 0.1882, 43.5523, 0.1882]) 

In [341]: io.imshow(img) 
Out[341]: <matplotlib.image.AxesImage at 0xce4d160>