On Jul 15, 2009, at 17:51, David Warde-Farley wrote: > Suppose I have an 8-bit integer 2-d array, X, and I want a 256x256 > matrix that tells me how many times a pixel value v was horizontally > (along second dimension) adjacent to a pixel value b
> Is there an efficient way to do such a thing with numpy operations? I > can't think of a way. I wrote the following cooccurrence matrix function as part of computing Haralick texture features. It's not exactly what you ask for, but close enough that it may give you some ideas. Vebjorn # From https://svnrepos.broad.mit.edu/CellProfiler/trunk/CellProfiler/pyCellProfiler/cellprofiler/cpmath/haralick.py import numpy as np import scipy.ndimage as scind def cooccurrence(quantized_image, labels, scale=3): """Calculates co-occurrence matrices for all the objects in the image. Return an array P of shape (nobjects, nlevels, nlevels) such that P[o, :, :] is the cooccurence matrix for object o. quantized_image -- a numpy array of integer type labels -- a numpy array of integer type scale -- an integer For each object O, the cooccurrence matrix is defined as follows. Given a row number I in the matrix, let A be the set of pixels in O with gray level I, excluding pixels in the rightmost S columns of the image. Let B be the set of pixels in O that are S pixels to the right of a pixel in A. Row I of the cooccurence matrix is the gray-level histogram of the pixels in B. """ nlevels = quantized_image.max() + 1 nobjects = labels.max() image_a = quantized_image[:, :-scale] image_b = quantized_image[:, scale:] labels_ab = labels[:, :-scale] equilabel = ((labels[:, :-scale] == labels[:, scale:]) & (labels[:,:-scale] > 0)) P, bins_P = np.histogramdd([labels_ab[equilabel]-1, image_a[equilabel], image_b[equilabel]], (nobjects, nlevels, nlevels)) pixel_count = fix(scind.sum(equilabel, labels[:,:-scale], np.arange(nobjects)+1)) pixel_count = np.tile(pixel_count[:,np.newaxis,np.newaxis], (1,nlevels,nlevels)) return (P.astype(float) / pixel_count.astype(float), nlevels) def fix(whatever_it_returned): """Convert a result from scipy.ndimage to a numpy array scipy.ndimage has the annoying habit of returning a single, bare value instead of an array if the indexes passed in are of length 1. For instance: scind.maximum(image, labels, [1]) returns a float but scind.maximum(image, labels, [1,2]) returns a list """ if getattr(whatever_it_returned,"__getitem__",False): return np.array(whatever_it_returned) else: return np.array([whatever_it_returned]) _______________________________________________ NumPy-Discussion mailing list NumPy-Discussion@scipy.org http://mail.scipy.org/mailman/listinfo/numpy-discussion
