用于GLCM计算和窗口大小的Mahotas库
我正在使用mahotas库对卫星图像(250×200像素)进行纹理分析(GLCM)。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。当代码完成后，我得到一个尺寸为6x8(而不是250×200)的图像，因为步骤大小已经设置为32。
所以，我的问题是:通过设置一个步长(为了避免在重叠区域中计算以及代码变得更快)，我是否可以以某种方式导出整个图像的GLCM结果，而不是它的子集(6x8维)？或者我别无选择，只能用正常的方式(不设置步长)遍历图像？
回答 1
Stack Overflow用户
发布于 2017-03-15 13:04:32
不能使用mahotas作为计算共现映射的函数，在这个库中不可用。从GLCM中提取纹理特征的另一种方法是利用skimage.feature.graycoprops (查看this thread获取细节)。
但是，如果您想坚持使用mahotas，您应该尝试使用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,))演示
对于下面运行的示例，我将win设置为19，它对应于形状(39, 39)的窗口。我考虑了两种不同的距离。注意，mht.haralick为四个方向生成了13个GLCM特性。总之，这将为每个像素生成一个104维特征向量。当应用于从陆地卫星图像中裁剪(250, 200)像素时，特征提取将在大约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>
https://stackoverflow.com/questions/42624729
复制相似问题
腾讯云开发者
