Schiebefenster in Python für GLCM Berechnung

Ich versuche, Texturanalyse in einem Satellitenbild mit GLCM-Algorithmus zu machen. Die Scikit-Image-Dokumentation ist dazu sehr hilfreich, aber für die GLCM-Berechnung benötigen wir eine Fenstergröße, die über das Bild schleift. Das ist zu langsam in Python. Ich fand viele Beiträge auf Stackoverflow über Schiebefenster, aber die Berechnung dauert ewig. Ich habe ein Beispiel unten gezeigt, es funktioniert aber für immer. Ich denke, das muss eine naive Art sein, es zu tun

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 

Ich kam mit skimage.util.view_as_windows Methode, die eine gute Lösung für mich sein könnte. Mein Problem ist, dass, wenn ich versuche, die GLCM zu berechnen, bekomme ich einen Fehler, der sagt:

ValueError: Das Parameterbild muss ein zweidimensionales Array sein

Dies liegt daran, dass das Ergebnis des GLCM-Bildes 4d Dimensionen hat und scikit-image view_as_windows Methode nur 2d Arrays akzeptiert. Hier ist mein Versuch

 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') 

Hat jemand eine Idee, wie ich die GLCM mit skimage.util.view_as_windows Methode berechnen kann?

One Solution collect form web for “Schiebefenster in Python für GLCM Berechnung”

Die Feature-Extraktion, die Sie durchführen möchten, ist eine computerintensive Aufgabe. Ich habe Ihre Methode beschleunigt, indem ich die Co-Vorkommen-Karte nur einmal für das gesamte Bild berechnet habe, anstatt die Co-Vorkommen-Karte immer wieder überlappende Positionen des Schiebefensters zu berechnen.

Die Co-Vorkommenskarte ist ein Stapel von Bildern gleicher Größe wie das Originalbild, bei dem – für jedes Pixel – Intensitätsniveau durch ganzzahlige Zahlen ersetzt wird, die das Co-Auftreten von zwei Intensitäten, nämlich Ii an diesem Pixel und Ij , codieren Bei einem Offset-Pixel. Die Co-Vorkommenskarte hat so viele Schichten, wie wir Offsets (dh alle möglichen Distanzwinkelpaare) betrachteten. Indem Sie die Co-Vorkommenskarte beibehalten müssen, müssen Sie die GLCM nicht an jeder Position des Schiebefensters vom Kratzer berechnen, da Sie die zuvor berechneten Co-Auftreten-Karten wiederverwenden können, um die Nachbarschaftsmatrizen (die GLCMs) für jede Distanz zu erhalten -angle-Paar Dieser Ansatz bietet Ihnen einen signifikanten Geschwindigkeitsgewinn.

Die Lösung, auf die ich kam, stützt sich auf die folgenden Funktionen:

 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 

DEMO

Die folgenden Ergebnisse entsprechen einem (250, 200) Pixel-Erntegut aus einem Landsat-Bild. Ich habe zwei Distanzen, vier Winkel und zwei GLCM-Eigenschaften berücksichtigt. Dies führt zu einem 16-dimensionalen Merkmalsvektor für jedes Pixel. Beachten Sie, dass das Schiebefenster quadriert ist und seine Seite 2*win + 1 Pixel (in diesem Test wurde ein Wert von win = 19 verwendet). Dieser Probelauf dauerte etwa 6 Minuten, was ziemlich kürzer ist als "für immer" 😉

 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> 

Satellitenbild

Python ist die beste Programmiersprache der Welt.