Реализация функции текстуры GLCM с помощью scikit-image и Python

Question

Я пытаюсь реализовать изображение текстуры как described in this tutorial с использованием Python и skimage.

Проблема заключается в перемещении окна 7x7 над большим растром и замене центра каждого пикселя на вычисленную текстуру из окна 7x7. Мне удается сделать это с помощью приведенного ниже кода, но я не вижу другого способа, кроме как прокручивать каждый отдельный пиксель, который очень медленный.

Один пакет программного обеспечения делает это через несколько секунд, поэтому должен быть какой-то другой способ ... есть?

Вот код, который работает, но очень медленно ...

import matplotlib.pyplot as plt 
import gdal, gdalconst 
import numpy as np 
from skimage.feature import greycomatrix, greycoprops 

filename = "//mnt//glaciology//RS2_20140101.jpg" 
outfilename = "//home//max//Documents//GLCM_contrast.tif" 
sarfile = gdal.Open(filename, gdalconst.GA_ReadOnly) 

sarraster = sarfile.ReadAsArray() 
#sarraster is satellite image, testraster will receive texture 
testraster = np.copy(sarraster) 
testraster[:] = 0 

for i in range(testraster.shape[0]): 
    print i, 
    for j in range(testraster.shape[1]): 

     #windows needs to fit completely in image 
     if i <3 or j <3: 
      continue 
     if i > (testraster.shape[0] - 4) or j > (testraster.shape[0] - 4): 
      continue 

     #Calculate GLCM on a 7x7 window 
     glcm_window = sarraster[i-3: i+4, j-3 : j+4] 
     glcm = greycomatrix(glcm_window, [1], [0], symmetric = True, normed = True) 

     #Calculate contrast and replace center pixel 
     contrast = greycoprops(glcm, 'contrast') 
     testraster[i,j]= contrast 

sarplot = plt.imshow(testraster, cmap = 'gray') 

Результаты:

Answer 1

У меня была такая же проблема, различные данные. Вот скрипт, я написал, что использует параллельную обработку и оконную подход скользящую:

import gdal, osr 
import numpy as np 
from scipy.interpolate import RectBivariateSpline 
from numpy.lib.stride_tricks import as_strided as ast 
import dask.array as da 
from joblib import Parallel, delayed, cpu_count 
import os 
from skimage.feature import greycomatrix, greycoprops 

def im_resize(im,Nx,Ny): 
    ''' 
    resize array by bivariate spline interpolation 
    ''' 
    ny, nx = np.shape(im) 
    xx = np.linspace(0,nx,Nx) 
    yy = np.linspace(0,ny,Ny) 

    try: 
     im = da.from_array(im, chunks=1000) #dask implementation 
    except: 
     pass 

    newKernel = RectBivariateSpline(np.r_[:ny],np.r_[:nx],im) 
    return newKernel(yy,xx) 

def p_me(Z, win): 
    ''' 
    loop to calculate greycoprops 
    ''' 
    try: 
     glcm = greycomatrix(Z, [5], [0], 256, symmetric=True, normed=True) 
     cont = greycoprops(glcm, 'contrast') 
     diss = greycoprops(glcm, 'dissimilarity') 
     homo = greycoprops(glcm, 'homogeneity') 
     eng = greycoprops(glcm, 'energy') 
     corr = greycoprops(glcm, 'correlation') 
     ASM = greycoprops(glcm, 'ASM') 
     return (cont, diss, homo, eng, corr, ASM) 
    except: 
     return (0,0,0,0,0,0) 


def read_raster(in_raster): 
    in_raster=in_raster 
    ds = gdal.Open(in_raster) 
    data = ds.GetRasterBand(1).ReadAsArray() 
    data[data<=0] = np.nan 
    gt = ds.GetGeoTransform() 
    xres = gt[1] 
    yres = gt[5] 

    # get the edge coordinates and add half the resolution 
    # to go to center coordinates 
    xmin = gt[0] + xres * 0.5 
    xmax = gt[0] + (xres * ds.RasterXSize) - xres * 0.5 
    ymin = gt[3] + (yres * ds.RasterYSize) + yres * 0.5 
    ymax = gt[3] - yres * 0.5 
    del ds 
    # create a grid of xy coordinates in the original projection 
    xx, yy = np.mgrid[xmin:xmax+xres:xres, ymax+yres:ymin:yres] 
    return data, xx, yy, gt 

def norm_shape(shap): 
    ''' 
    Normalize numpy array shapes so they're always expressed as a tuple, 
    even for one-dimensional shapes. 
    ''' 
    try: 
     i = int(shap) 
     return (i,) 
    except TypeError: 
     # shape was not a number 
     pass 

    try: 
     t = tuple(shap) 
     return t 
    except TypeError: 
     # shape was not iterable 
     pass 

    raise TypeError('shape must be an int, or a tuple of ints') 

def sliding_window(a, ws, ss = None, flatten = True): 
    ''' 
    Source: http://www.johnvinyard.com/blog/?p=268#more-268 
    Parameters: 
     a - an n-dimensional numpy array 
     ws - an int (a is 1D) or tuple (a is 2D or greater) representing the size 
      of each dimension of the window 
     ss - an int (a is 1D) or tuple (a is 2D or greater) representing the 
      amount to slide the window in each dimension. If not specified, it 
      defaults to ws. 
     flatten - if True, all slices are flattened, otherwise, there is an 
        extra dimension for each dimension of the input. 

    Returns 
     an array containing each n-dimensional window from a 
    '''  
    if None is ss: 
     # ss was not provided. the windows will not overlap in any direction. 
     ss = ws 
    ws = norm_shape(ws) 
    ss = norm_shape(ss) 
    # convert ws, ss, and a.shape to numpy arrays 
    ws = np.array(ws) 
    ss = np.array(ss) 
    shap = np.array(a.shape) 
    # ensure that ws, ss, and a.shape all have the same number of dimensions 
    ls = [len(shap),len(ws),len(ss)] 
    if 1 != len(set(ls)): 
     raise ValueError(\ 
     'a.shape, ws and ss must all have the same length. They were %s' % str(ls)) 

    # ensure that ws is smaller than a in every dimension 
    if np.any(ws > shap): 
     raise ValueError(\ 
     'ws cannot be larger than a in any dimension.\ 
    a.shape was %s and ws was %s' % (str(a.shape),str(ws))) 

    # how many slices will there be in each dimension? 
    newshape = norm_shape(((shap - ws) // ss) + 1) 


    # the shape of the strided array will be the number of slices in each dimension 
    # plus the shape of the window (tuple addition) 
    newshape += norm_shape(ws) 


    # the strides tuple will be the array's strides multiplied by step size, plus 
    # the array's strides (tuple addition) 
    newstrides = norm_shape(np.array(a.strides) * ss) + a.strides 
    a = ast(a,shape = newshape,strides = newstrides) 
    if not flatten: 
     return a 
    # Collapse strided so that it has one more dimension than the window. I.e., 
    # the new array is a flat list of slices. 
    meat = len(ws) if ws.shape else 0 
    firstdim = (np.product(newshape[:-meat]),) if ws.shape else() 
    dim = firstdim + (newshape[-meat:]) 
    # remove any dimensions with size 1 
    dim = filter(lambda i : i != 1,dim) 

    return a.reshape(dim), newshape 

def CreateRaster(xx,yy,std,gt,proj,driverName,outFile): 
    ''' 
    Exports data to GTiff Raster 
    ''' 
    std = np.squeeze(std) 
    std[np.isinf(std)] = -99 
    driver = gdal.GetDriverByName(driverName) 
    rows,cols = np.shape(std) 
    ds = driver.Create(outFile, cols, rows, 1, gdal.GDT_Float32)  
    if proj is not None: 
     ds.SetProjection(proj.ExportToWkt()) 
    ds.SetGeoTransform(gt) 
    ss_band = ds.GetRasterBand(1) 
    ss_band.WriteArray(std) 
    ss_band.SetNoDataValue(-99) 
    ss_band.FlushCache() 
    ss_band.ComputeStatistics(False) 
    del ds 


#Stuff to change 

if __name__ == '__main__': 
    win_sizes = [7] 
    for win_size in win_sizes[:]: 
     in_raster = #Path to input raster 
     win = win_size 
     meter = str(win/4) 

     #Define output file names 
     contFile = 
     dissFile = 
     homoFile = 
     energyFile = 
     corrFile = 
     ASMFile = 



     merge, xx, yy, gt = read_raster(in_raster) 

     merge[np.isnan(merge)] = 0 

     Z,ind = sliding_window(merge,(win,win),(win,win)) 

     Ny, Nx = np.shape(merge) 

     w = Parallel(n_jobs = cpu_count(), verbose=0)(delayed(p_me)(Z[k]) for k in xrange(len(Z))) 

     cont = [a[0] for a in w] 
     diss = [a[1] for a in w] 
     homo = [a[2] for a in w] 
     eng = [a[3] for a in w] 
     corr = [a[4] for a in w] 
     ASM = [a[5] for a in w] 


     #Reshape to match number of windows 
     plt_cont = np.reshape(cont , (ind[0], ind[1])) 
     plt_diss = np.reshape(diss , (ind[0], ind[1])) 
     plt_homo = np.reshape(homo , (ind[0], ind[1])) 
     plt_eng = np.reshape(eng , (ind[0], ind[1])) 
     plt_corr = np.reshape(corr , (ind[0], ind[1])) 
     plt_ASM = np.reshape(ASM , (ind[0], ind[1])) 
     del cont, diss, homo, eng, corr, ASM 

     #Resize Images to receive texture and define filenames 
     contrast = im_resize(plt_cont,Nx,Ny) 
     contrast[merge==0]=np.nan 
     dissimilarity = im_resize(plt_diss,Nx,Ny) 
     dissimilarity[merge==0]=np.nan  
     homogeneity = im_resize(plt_homo,Nx,Ny) 
     homogeneity[merge==0]=np.nan 
     energy = im_resize(plt_eng,Nx,Ny) 
     energy[merge==0]=np.nan 
     correlation = im_resize(plt_corr,Nx,Ny) 
     correlation[merge==0]=np.nan 
     ASM = im_resize(plt_ASM,Nx,Ny) 
     ASM[merge==0]=np.nan 
     del plt_cont, plt_diss, plt_homo, plt_eng, plt_corr, plt_ASM 


     del w,Z,ind,Ny,Nx 

     driverName= 'GTiff'  
     epsg_code=26949 
     proj = osr.SpatialReference() 
     proj.ImportFromEPSG(epsg_code) 

     CreateRaster(xx, yy, contrast, gt, proj,driverName,contFile) 
     CreateRaster(xx, yy, dissimilarity, gt, proj,driverName,dissFile) 
     CreateRaster(xx, yy, homogeneity, gt, proj,driverName,homoFile) 
     CreateRaster(xx, yy, energy, gt, proj,driverName,energyFile) 
     CreateRaster(xx, yy, correlation, gt, proj,driverName,corrFile) 
     CreateRaster(xx, yy, ASM, gt, proj,driverName,ASMFile) 

     del contrast, merge, xx, yy, gt, meter, dissimilarity, homogeneity, energy, correlation, ASM 

Этот скрипт вычисляет свойства КРНБА для определенного размера окна, без перекрытия между соседними окнами.