在Python中滑动window进行GLCM计算
Sliding window in Python for GLCM calculation
我正在尝试使用 GLCM 算法在卫星图像中进行纹理分析。 scikit-image 文档对此非常有帮助，但对于 GLCM 计算，我们需要 window 大小的图像循环。这在 Python 中太慢了。我在 Whosebug 上发现了很多关于滑动 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
我遇到了 skimage.util.view_as_windows 方法，这对我来说可能是一个很好的解决方案。我的问题是，当我尝试计算 GLCM 时，我收到一条错误消息：
ValueError: The parameter image must be a 2-dimensional array
这是因为 GLCM 图像的结果具有 4d 维度，而 scikit-image 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 吗？
您尝试执行的特征提取是一项计算机密集型任务。我通过对整个图像只计算一次 共现图 来加快你的方法，而不是在滑动 [=36] 的重叠位置上一遍又一遍地计算共现图=].
共现图是一堆与原始图像大小相同的图像，其中 - 对于每个像素 - 强度级别被编码两个强度共现的整数代替，即 Ii 在该像素处，Ij 在偏移像素处。共现图的层数与我们考虑的偏移量（即所有可能的距离-角度对）一样多。通过保留共现图，你不需要从头开始计算滑动window的每个位置的GLCM，因为你可以重复使用之前计算的共现图来获得邻接矩阵（GLCMs ) 对于每个距离-角度对。这种方法可以显着提高速度。
我想出的解决方案依赖于以下功能：
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
演示版
以下结果对应于 Landsat 图像的 (250, 200) 像素裁剪。我考虑了两个距离、四个角度和两个 GLCM 属性。这导致每个像素的 16 维特征向量。请注意，滑动 window 是正方形，它的边是 2*win + 1 像素（在此测试中使用了 win = 19 的值）。此示例 运行 花费了大约 6 分钟，比 "forever" ;-)
短得多
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>
我正在尝试使用 GLCM 算法在卫星图像中进行纹理分析。 scikit-image 文档对此非常有帮助，但对于 GLCM 计算，我们需要 window 大小的图像循环。这在 Python 中太慢了。我在 Whosebug 上发现了很多关于滑动 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
我遇到了 skimage.util.view_as_windows 方法，这对我来说可能是一个很好的解决方案。我的问题是，当我尝试计算 GLCM 时，我收到一条错误消息：
ValueError: The parameter
imagemust be a 2-dimensional array
这是因为 GLCM 图像的结果具有 4d 维度，而 scikit-image 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 吗？
您尝试执行的特征提取是一项计算机密集型任务。我通过对整个图像只计算一次 共现图 来加快你的方法，而不是在滑动 [=36] 的重叠位置上一遍又一遍地计算共现图=].
共现图是一堆与原始图像大小相同的图像，其中 - 对于每个像素 - 强度级别被编码两个强度共现的整数代替，即 Ii 在该像素处，Ij 在偏移像素处。共现图的层数与我们考虑的偏移量（即所有可能的距离-角度对）一样多。通过保留共现图，你不需要从头开始计算滑动window的每个位置的GLCM，因为你可以重复使用之前计算的共现图来获得邻接矩阵（GLCMs ) 对于每个距离-角度对。这种方法可以显着提高速度。
我想出的解决方案依赖于以下功能：
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
演示版
以下结果对应于 Landsat 图像的 (250, 200) 像素裁剪。我考虑了两个距离、四个角度和两个 GLCM 属性。这导致每个像素的 16 维特征向量。请注意，滑动 window 是正方形，它的边是 2*win + 1 像素（在此测试中使用了 win = 19 的值）。此示例 运行 花费了大约 6 分钟，比 "forever" ;-)
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>