[Numpy-discussion] Neighbour-frequency matrix
Wed Jul 15 17:40:54 CDT 2009
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.
# 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
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,
(nobjects, nlevels, nlevels))
pixel_count = fix(scind.sum(equilabel, labels[:,:-scale],
pixel_count = np.tile(pixel_count[:,np.newaxis,np.newaxis],
return (P.astype(float) / pixel_count.astype(float), nlevels)
"""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.
scind.maximum(image, labels, ) returns a float
scind.maximum(image, labels, [1,2]) returns a list
More information about the NumPy-Discussion