Fenêtre coulissante en Python pour le calcul GLCM

Question

J'essaie de faire une analyse de texture dans une imagerie satellite en utilisant l'algorithme GLCM. La documentation de scikit-image est très utile à ce sujet mais pour le calcul de GLCM nous avons besoin d'une taille de fenêtre en boucle sur l'image. C'est trop lent en Python. J'ai trouvé de nombreux posts sur stackoverflow à propos des fenêtres coulissantes mais le calcul prend pour toujours. J'ai un exemple ci-dessous, cela fonctionne mais prend une éternité. Je suppose que cela doit être une façon naïve de le faire

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 

je suis tombé avec la méthode skimage.util.view_as_windows qui pourrait être une bonne solution pour moi. Mon problème est que, lorsque je tente de calculer le GLCM je reçois une erreur qui dit:

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

En effet, le résultat de la méthode image GLCM a des dimensions 4d et scikit image view_as_windows accepte des tableaux 2D seulement. Voici ma tentative

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') 

Est-ce que quelqu'un a une idée sur la façon dont je peux calculer la GLCM en utilisant la méthode skimage.util.view_as_windows?

Answer 1

L'extraction de fonctions que vous essayez d'effectuer est une tâche exigeante en informatique. J'ai accéléré votre méthode en calculant la carte de cooccurrence une seule fois pour l'ensemble de l'image, plutôt que de calculer la carte de cooccurrence à plusieurs reprises sur les positions de chevauchement de la fenêtre glissante. La carte de cooccurrence est une pile d'images de la même taille que l'image originale, dans laquelle - pour chaque pixel - les niveaux d'intensité sont remplacés par des nombres entiers qui codent la cooccurrence de deux intensités, à savoir Ii à ce pixel et Ij à un pixel décalé. La carte de cooccurrence a autant de couches que nous avons considérées de décalages (c'est-à-dire toutes les paires distance-angle possibles). En conservant la carte de cooccurrence, vous n'avez pas besoin de calculer le GLCM à chaque position de la fenêtre glissante, car vous pouvez réutiliser les cartes de cooccurrence précédemment calculées pour obtenir les matrices d'adjacence (les GLCM) pour chaque distance paire -angle. Cette approche vous fournit un gain de vitesse significatif.

La solution que je suis venu avec se fie aux fonctions suivantes:

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

Les résultats suivants correspondent à une culture (250, 200) pixels d'une image Landsat. J'ai considéré deux distances, quatre angles et deux propriétés GLCM. Cela se traduit par un vecteur de caractéristiques en 16 dimensions pour chaque pixel. Notez que la fenêtre glissante est au carré et que son côté est 2*win + 1 pixels (dans ce test, une valeur de win = 19 a été utilisée). Cet échantillon analysé a environ 6 minutes, ce qui est assez courte que « jamais » ;-)

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>