Neighboring gray-level dependence matrix (NGLDM) in MATLAB

Question

I would like to calculate a couple of texture features (namely: small/ large number emphasis, number non-uniformity, second moment and entropy). Those can be computed from Neighboring gray-level dependence matrix. I'm struggling with understanding/implementation of this. There is very little info on this method (publicly available).

According to this paper:

This matrix takes the form of a two-dimensional array Q, where Q(i,j) can be considered as frequency counts of grayness variation of a processed image. It has a similar meaning as histogram of an image. This array is N_g×N_r where N_g is the number of possible gray levels and N_r is the number of possible neighbours to a pixel in an image.

If the image function f(i,j) is discrete, then it is easy to computer the Q matrix (for positive integer d, a) by counting the number of times the difference between each element in f(i,j) and its neighbours is equal or less than a at a certain distance d.

Here is the example from the same paper (d = 1, a = 0):

Input (image) matrix and output matrix Q:

I've been looking at this example for hours now and still can't figure out how they got that Q matrix. Anyone?

_{The method was originally created by C. Sun and W. Wee and was described in a paper called: "Neighboring gray level dependence matrix for texture classification" to which I got access, but can't download (after pressing download the page reloads and that's it).}

Answer 1

In the example that you have provided, d=1 and a=0. When d=1, we consider pixels in an 8-pixel neighbourhood. When a=0, this means that we look for pixels that have the same value as the centre of the neighbourhood.

The basic algorithm is the following:

1. Initialize your NGLDM matrix to all zeroes. The total number of rows corresponds to the total number of possible intensities / values in your image. The total number of columns corresponds to how many pixels are in your neighbourhood plus 1. As such for d=1, we have an 8-pixel neighbourhood and so 8 + 1 = 9. Because there are 4 possible intensities (0,1,2,3), we thus have a 4 x 9 matrix. Let's call this matrix M.
2. For each pixel in your matrix, take note of this pixel. This goes in the Ng row.
3. Write out how many valid neighbours there are that surround this pixel.
4. Count how many times you see the neighbouring pixels matching that pixel in Step #1. This is your Nr column.
5. Once you figure out the numbers in Step #1 and Step #2, increment this location by 1.

---

Here's a slight gotcha: They ignore the border locations. As such, you don't do this procedure for the first row, last row, first column or last column. My guess is that they want to be sure that you have an 8-pixel neighbourhood all the time. This is also dictated by the distance d=1. You must be able to grab every valid pixel given a centre location at d=1. If d=2, then you would have to make sure that every pixel in the centre of the neighbourhood has a 25 pixel neighbourhood and so on.

Let's start from the second row, second column location of this matrix. Let's go through the steps:

1. Ng = 1 as the location is 1.
2. Valid neighbours - Starting from the top left pixel in this neighbourhood, and scanning left to right and omitting the centre, we have: 1, 1, 2, 0, 1, 0, 2, 2.
3. How many values are equal to 1? Three times. Therefore Nr = 3
4. M(Ng,Nr) += 1. Access row Ng = 1, and access row Nr = 3, and increment this spot by 1.

---

Want to know how I figured out they don't count the borders? Let's do the bottom left pixel. That location is 0, so Ng = 0. If you repeat the algorithm that I just said, you would expect Ng = 0, Nr = 1, and so you would expect at least one entry in that location in your matrix... but you don't! If you do similar checks around the border of the image, you'll see that entries that are supposed to be there... aren't. Take a look at the third row, fifth column. You would think that Ng = 1 and Nr = 1, but we don't see that in the matrix.

---

One more example. Why is M(Ng,Nr) = 4, Ng = 2, Nr = 4? Well, take a look at every pixel that has a 2 in it. The only valid locations where we can capture an 8 pixel neighbourhood successfully are the row=2, col=4, row=3, col=3, row=3, col=4, row=4, col=3, and row=4, col=4. By applying the same algorithm that we have seen, you'll see that for each of those locations, Nr = 4. As such, we see this combination of Ng = 2, Nr = 4 four times, and that's why the location is set to 4. However, in row=3, col=4, this actually is Nr = 5, as there are five 2s in that neighbourhood at that centre. That's why you see Ng = 2, Nr = 5, M(Ng,Nr) = 1.

As an example, let's do one of the locations. Let's do the 2 smack dab in the middle of the matrix (row=3, col=3):

1. Ng = 2
2. What are the valid neighbouring pixels? 1, 1, 2, 0, 2, 3, 2, 2 (omit the centre)
3. Count how many pixels equal to 2. There are four of them, so Nr = 4
4. M(Ng,Nr) += 1. Take Ng = 2, Nr = 4 and increment this spot by 1.

If you do this with the other valid locations that have 2, you'll see that Nr = 4 each time with the exception of the third row and fourth column, where Nr = 5.

---

So how would we implement this in MATLAB? What you can do is use im2col to transform each valid neighbourhood into columns. What I'm also going to do is extract the centre of each neighbourhood. This is actually the middle row of the matrix. We will then figure out how many pixels for each neighbourhood equal the centre, sum them up, and this will determine our Nr values. The Ng values will be the middle row values themselves. Once we do this, we can compute a histogram based on these values just like how the algorithm is doing to get our matrix. In other words, try doing this:

% // Your example
A = [1 1 2 3 1; 0 1 1 2 2; 0 0 2 2 1; 3 3 2 2 1; 0 0 2 0 1];
B = im2col(A, [3 3]); %//Convert neighbourhoods to columns - 3 x 3 means d = 1
C = bsxfun(@eq, B, B(5,:)); %//Figure out a logical matrix where each column tells
                            %//you how many elements equals the one in each centre
D = sum(C, 1) - 1; %// Must subtract by 1 to discount centre pixel
Ng = B(5,:).' + 1; % // We must make this into a column vector, and we also must
                   % // offset by 1 as MATLAB starts indexing by 1.
                   %// Column vector is for accumarray input
Nr = D.' + 1; %// Do the same for Nr.  We could have simply left out the + 1 here and
              %// took out the subtraction of -1 for D, but I want to explicitly show
              %// the steps
Q = accumarray([Ng Nr], 1, [4 9]); %// 4 unique intensities, 9 possible locations (0-8)

... and here is our matrix:

Q =

 0     0     1     0     0     0     0     0     0
 0     0     1     1     0     0     0     0     0
 0     0     0     0     4     1     0     0     0
 0     1     0     0     0     0     0     0     0

If you check this, you'll see this matches with Q.

---

Bonus

If you want to be able to accommodate for the algorithm in general, where you specify d and a, we can simply follow the guidelines of your text. For each neighbourhood, you find the difference between the centre pixel and all of the other pixels. You count how many pixels are <= a for any positive integer d. Note that this will create a 2*d + 1 x 2*d + 1 neighbourhood we need to examine. We can also make this into a function. Without further ado:

%// Set A up yourself, then use a and d as inputs
%// Precondition - a and d are both integers. a can be 0 and d is positive!
function [Q] = calculateGrayDepMatrix(A, a, d)
    neigh = 2*d + 1; % //Calculate rows/columns of neighbourhood
    numTotalNeigh = neigh*neigh; % //Calculate total number of pixels in neighbourhood
    middleRow = ceil(numTotalNeigh / 2); %// Figure out which index the middle row is
    B = im2col(A, [neigh neigh]);  %// Make into columns
    Cdiff = abs(bsxfun(@minus, B, B(middleRow,:))); %// For each neighbourhood, subtract with its centre
    C = Cdiff <= a; %// For each neighbourhood, figure out which differences are <= a
    D = sum(C, 1) - 1; % //For each neighbourhood, add them up
    Ng = B(middleRow,:).' + 1; % // Determine Ng and Nr, and find Q
    Nr = D.' + 1; 
    Q = accumarray([Ng Nr], 1, [max(Ng) numTotalNeigh]); 
end

We can recreate the scenario we showed above with the example matrix by:

A = [1 1 2 3 1; 0 1 1 2 2; 0 0 2 2 1; 3 3 2 2 1; 0 0 2 0 1];
Q = calculateGrayDepMatrix(A, 0, 1);

Q is thus:

Q =

 0     0     1     0     0     0     0     0     0
 0     0     1     1     0     0     0     0     0
 0     0     0     0     4     1     0     0     0
 0     1     0     0     0     0     0     0     0

Hope this helps!