Скользящее окно в Python для расчета GLCM

Question

Я пытаюсь сделать анализ текстуры на спутниковых снимках, используя алгоритм GLCM. Документация scikit-image очень полезна для этого, но для вычисления GLCM нам нужен размер окна, зацикливающийся на изображении. Это слишком медленно в Python. Я нашел много сообщений на stackru о скользящих окнах, но вычисления требуют вечно. У меня есть пример, показанный ниже, он работает, но работает вечно. Я предполагаю, что это должно быть наивным способом сделать это

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 

Я сталкивался с skimage.util.view_as_windows метод, который может быть хорошим решением для меня. Моя проблема в том, что, когда я пытаюсь вычислить GLCM, я получаю сообщение об ошибке:

ValueError: параметр image должен быть двумерным массивом

Это связано с тем, что результат изображения 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')

Кто-нибудь есть идеи о том, как я могу рассчитать GLCM с помощью skimage.util.view_as_windows метод?

Answer 1

Извлечение функции, которую вы пытаетесь выполнить, - задача, требующая большого количества компьютеров. Я ускорил ваш метод, вычисляя карту совместного появления только один раз для всего изображения, а не вычисляя карту совместного появления снова и снова на перекрывающихся позициях скользящего окна.

Карта совместного вхождения представляет собой стек изображений того же размера, что и исходное изображение, в котором - для каждого пикселя - уровни интенсивности заменяются целыми числами, которые кодируют совместное вхождение двух интенсивностей, а именно Ii в этом пикселе и Ij на смещенном пикселе. Карта совместного вхождения имеет столько слоев, сколько мы рассмотрели смещения (т. Е. Все возможные пары расстояние-угол). Сохраняя карту совместного появления, вам не нужно вычислять GLCM в каждой позиции скользящего окна с нуля, так как вы можете повторно использовать ранее вычисленные карты совместного появления, чтобы получить матрицы смежности (GLCM) для каждого расстояния пара Такой подход дает вам значительное увеличение скорости.

Решение, которое я придумал, опирается на функции ниже:

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

Следующие результаты соответствуют (250, 200) пиксели обрезаются с изображения Landsat. Я рассмотрел два расстояния, четыре угла и два свойства GLCM. В результате получается 16-мерный вектор признаков для каждого пикселя. Обратите внимание, что скользящее окно имеет квадрат и его сторона 2*win + 1 пикселей (в этом тесте значение win = 19 использовался). Этот примерный прогон занял около 6 минут, что короче, чем "навсегда";-)

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>