[Numpy-discussion] Neighbour-frequency matrix

[Vebjorn Ljosa]

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