Python: получение GLCM непрямоугольной области

Хотя mahotas также является отличной библиотекой компьютерного зрения, для этого нет необходимости прекращать использовать skimage.

Что необходимо необходимо, как указал @Tonechas, так это установить для этих значений NaN целое число, поскольку np.nan имеет тип float, а для функции greycomatrix требуется массив целых чисел.

Самым простым вариантом было бы установить эти NaN на ноль, но, если у вас уже есть нулевые значения в ваших пикселях и вы не хотите их смешивать, вы можете выбрать любую другую константу. После этого все, что вам нужно сделать, это отфильтровать выбранное значение (опять же, как правило, ноль) из GLCM.

Чтобы понять, что это значит, давайте посмотрим, что skimage говорит нам о выводе greycomatrixфункция:

4D ndarray

[...] Значение P[i,j,d,theta] - это количество раз, когда уровень серого j встречается на расстоянии d и под углом тета от уровня серого i. Если normed имеет значение False, выходные данные имеют тип uint32, в противном случае — float64. Размеры: уровни x уровни x количество расстояний x количество углов.

Другими словами, первые два измерения массива определяют матрицу, которая сообщает нам, сколько раз два разных значения находятся на определенном расстоянии друг от друга. Обратите внимание, что GLCM не сохраняет форму входного массива. Эти строки и столбцы говорят нам, как соотносятся значения.

Зная это, легко отфильтровать значения за пределами нашего ROI (представьте, что мы установили эти значения NaN равными нулю):

glcm = greycomatrix(img, [1], [0])  # Calculate the GLCM "one pixel to the right"
filt_glcm = glcm[1:, 1:, :, :]           # Filter out the first row and column

Теперь вы можете легко вычислить свойства Haralick вашего отфильтрованного GLCM. Например:

greycoprops(filt_glcm, prop='contrast')

Проблема в том, что вам нужно передать целочисленный массив в greycomatrix, но np.nan имеет тип float (взгляните на этот поток для подробностей). В результате вы не можете кодировать пиксели за пределами ROI как NaN.

Приблизительным обходным путем для работы с непрямоугольными ROI будет установка пикселей за пределами ROI на 0 и использование функции haralick из mahotas. Эта функция возвращает 13 признаков Харалика, извлеченных из четырех разных GLCM, соответствующих четырем двумерным ориентациям и конкретному значению параметра расстояния.

Из документации:

ignore_zeros можно использовать, чтобы функция игнорировала любые пиксели с нулевым значением (в качестве фона).

Таким образом, вам нужно замаскировать те пиксели, которые выходят за пределы ROI, и установить ignore_zeros в True в вызове haralick.

ДЕМО

Для начала давайте сгенерируем некоторые фиктивные данные:

In [213]: import numpy as np

In [214]: shape = (3, 4)

In [215]: levels = 8

In [216]: np.random.seed(2017)

In [217]: x = np.random.randint(0, levels, size=shape)

In [218]: x
Out[218]: 
array([[3, 1, 6, 5],
       [2, 0, 2, 2],
       [3, 7, 7, 7]])

Затем мы должны удалить все нули из изображения, так как в этом подходе нулевой уровень интенсивности зарезервирован для пикселей за пределами области интереса. Стоит отметить, что объединение интенсивностей 0 и 1 в одну интенсивность 1 вносит неточности в результаты.

In [219]: x[x == 0] = 1

In [220]: x
Out[220]: 
array([[3, 1, 6, 5],
       [2, 1, 2, 2],
       [3, 7, 7, 7]])

Следующий шаг состоит в определении маски для пикселей за пределами области интереса (в этом игрушечном примере — четырех углах изображения) и установке для этих пикселей значения 0.

In [221]: non_roi = np.zeros(shape=shape, dtype=np.bool)

In [222]: non_roi[np.ix_([0, -1], [0, -1])] = True

In [223]: non_roi
Out[223]: 
array([[ True, False, False,  True],
       [False, False, False, False],
       [ True, False, False,  True]], dtype=bool)

In [224]: x[non_roi] = 0

In [225]: x
Out[225]: 
array([[0, 1, 6, 0],
       [2, 1, 2, 2],
       [0, 7, 7, 0]])

Теперь мы можем выполнить извлечение признаков из GLCM непрямоугольной области интереса:

In [226]: import mahotas.features.texture as mht

In [227]: features = mht.haralick(x, ignore_zeros=True)

In [228]: features.size
Out[228]: 52

In [229]: features.ravel()
Out[229]: array([ 0.18      ,  5.4       ,  0.5254833 , ...,  0.81127812,
       -0.68810414,  0.96300727])

Может быть полезно проверить, как выглядят матрицы совпадений. Например, GLCM «пиксель вправо» будет:

In [230]: mht.cooccurence(x, 0)
Out[230]: 
array([[0, 1, 0, ..., 0, 1, 2],
       [1, 0, 2, ..., 0, 1, 0],
       [0, 2, 2, ..., 0, 0, 0],
       ..., 
       [0, 0, 0, ..., 0, 0, 0],
       [1, 1, 0, ..., 0, 0, 0],
       [2, 0, 0, ..., 0, 0, 2]])