[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.


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

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
     scind.maximum(image, labels, [1,2]) returns a list
     if getattr(whatever_it_returned,"__getitem__",False):
         return np.array(whatever_it_returned)
         return np.array([whatever_it_returned])

More information about the NumPy-Discussion mailing list