Библиотека Mahotas для вычисления GLCM и размера окна

Question

Я использую библиотеку mahotas для анализа текстур (GLCM) на спутниковом изображении (250 x 200 пикселей). Расчеты GLCM выполняются в пределах размера окна. Таким образом, для двух соседних позиций скользящего окна нам нужно с нуля вычислить две матрицы совместного использования. Я прочитал, что я могу также установить размер шага, чтобы избежать вычисления GLCM в перекрывающихся областях. Я предоставил код ниже.

#Compute haralick features
def haralick_feature(image):
    haralick = mahotas.features.haralick(image, True)
    return haralick


img = 'SAR_image.tif'
win_s=32 #window size
step=32 #step size

rows = img.shape[0]
cols = img.shape[1]
array = np.zeros((rows,cols), dtype= object)
harList = []

for i in range(0, rows-win_s-1, step):
        print 'Row number: ', r
    for j in range(0, cols-win_s-1, step):
        harList.append(haralick_feature(image))

harImages = np.array(harList)     
harImages_mean = harImages.mean(axis=1)

Для приведенного выше кода я установил размер окна и шага равным 32. Когда код завершится, я получу изображение с размерами 6 x 8 (вместо 250 x 200), это имеет смысл, так как размер шага был установлен равным 32,

Итак, мой вопрос: установив размер шага (чтобы избежать вычислений в перекрывающихся областях, а также кода становится быстрее), я могу каким-то образом получить результаты GLCM для всего изображения с размерами 250 x 200 вместо того, чтобы иметь его подмножество (6 х 8 размеров) Или у меня нет выбора, кроме как зацикливание изображения обычным способом (без установки размера шага)?

Answer 1

Вы не можете сделать это с помощью mahotas поскольку функция, которая вычисляет карту совместного использования, не доступна в этой библиотеке. Альтернативный подход к извлечению текстурных особенностей из GLCM состоит в использовании skimage.feature.graycoprops (проверьте эту тему для деталей).

Но если вы хотите придерживаться махот, вы должны попробовать использовать skimage.util.view_as_windows а не скользящее окно, поскольку это может ускорить сканирование по изображению. Пожалуйста, не забудьте прочитать предостережение об использовании памяти в конце документации. При использовании view_as_windows это доступный для вас подход, следующий код выполняет свою работу:

import numpy as np
from skimage import io, util
import mahotas.features.texture as mht

def haralick_features(img, win, d):
    win_sz = 2*win + 1
    window_shape = (win_sz, win_sz)
    arr = np.pad(img, win, mode='reflect')
    windows = util.view_as_windows(arr, window_shape)
    Nd = len(d)
    feats = np.zeros(shape=windows.shape[:2] + (Nd, 4, 13), dtype=np.float64)
    for m in xrange(windows.shape[0]):
        for n in xrange(windows.shape[1]):
            for i, di in enumerate(d):
                w = windows[m, n, :, :]
                feats[m, n, i, :, :] = mht.haralick(w, distance=di)
    return feats.reshape(feats.shape[:2] + (-1,))

---

DEMO

Для примера запуска ниже я установил win в 19 который соответствует окну формы (39, 39), Я рассмотрел два разных расстояния. Заметить, что mht.haralick дает 13 GLCM функций для четырех направлений. В целом это приводит к 104-мерному вектору признаков для каждого пикселя. Применительно к (250, 200) пиксели обрезаются с изображения Landsat, извлечение объекта заканчивается примерно через 7 минут.

In [171]: img = io.imread('landsat_crop.tif')

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

In [173]: win = 19

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

In [175]: %time feature_map = haralick_features(img, win, d)
Wall time: 7min 4s

In [176]: feature_map.shape
Out[176]: (250L, 200L, 104L)

In [177]: feature_map[0, 0, :]
Out[177]: 
array([  8.19278030e-03,   1.30863698e+01,   7.64234582e-01, ...,
         3.59561817e+00,  -1.35383606e-01,   8.32570045e-01])

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