# [Numpy-discussion] Neighbour-frequency matrix

Vebjorn Ljosa vebjorn@ljosa....
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.

Vebjorn

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

```