MATLAB 中的相邻灰度依赖矩阵 (NGLDM)

Question

我想计算几个纹理特征（即：小/大数字强调、数字非均匀性、二阶矩和熵）。这些可以从相邻灰度级依赖矩阵计算出来。我正在努力理解/实现这一点。关于这种方法的信息很少（公开可用）。

根据这篇论文：

该矩阵采用二维数组Q 的形式，其中Q(i,j)可以被视为已处理图像灰度变化的频率计数。它与图像的直方图具有相似的含义。该数组为 N _g ×N _r，其中 N _g是可能的灰度级数，N _r是图像中像素的可能邻居数。

如果图像函数F（I，J）是离散的，则很容易计算机的Q矩阵（对于正整数d，一）通过在计数的次数每个元素之间的差F（I，J）和其在特定距离d处，邻居等于或小于a。

这是同一篇论文中的示例（d = 1, a = 0）：

输入（图像）矩阵和输出矩阵Q：

我已经看了这个例子好几个小时了，但仍然无法弄清楚他们是如何得到Q矩阵的。任何人？

_{该方法最初由 C. Sun 和 W. Wee 创建，并在一篇名为“纹理分类的相邻灰度依赖矩阵”的论文中进行了描述，我可以访问，但无法下载（按下载后页面重新加载和就是这样）。}

Answer 1

在您提供的示例中，d=1和a=0. 当 时d=1，我们考虑 8 像素邻域中的像素。当 时a=0，这意味着我们寻找与邻域中心具有相同值的像素。

基本算法如下：

1. 将 NGLDM 矩阵初始化为全零。总行数对应于图像中可能的强度/值的总数。总列数对应于您的邻域中的像素数加 1。因此d=1，我们有一个 8 像素的邻域，因此 8 + 1 = 9。因为有 4 种可能的强度(0,1,2,3)，因此我们有一个 4 x 9矩阵。我们称这个矩阵为M。
2. 对于矩阵中的每个像素，记下这个像素。这是Ng排成一排的。
3. 写出多少有效的邻居也有围绕该像素。
4. 计算您在步骤 1 中看到与该像素匹配的相邻像素的次数。这是你的Nr专栏。
5. 计算出步骤 #1 和步骤 #2 中的数字后，将此位置增加 1。

这里有一个小问题： 他们忽略了边界位置。因此，您不会对第一行、最后一行、第一列或最后一列执行此过程。我的猜测是他们希望确保您始终拥有一个 8 像素的邻域。这也是由距离决定的d=1。您必须能够抓取给定中心位置的每个有效像素d=1。如果d=2，那么您必须确保邻域中心的每个像素都有一个 25 像素的邻域，依此类推。

让我们从这个矩阵的第二行第二列位置开始。让我们来看看这些步骤：

1. Ng = 1 因为位置是 1。
2. 有效邻-从顶部开始在这附近离开像素和扫描从左到右省略中心，我们有：1, 1, 2, 0, 1, 0, 2, 2。
3. 等于多少个值1？三次。所以Nr = 3
4. M(Ng,Nr) += 1. 访问 rowNg = 1和 access row Nr = 3，并将该位置加 1。

想知道我是怎么知道他们不计算边界的吗？让我们做左下角的像素。那个位置是0，所以Ng = 0。如果您重复我刚才所说的算法，您会期望Ng = 0, Nr = 1，因此您会期望矩阵中该位置至少有一个条目……但您没有！如果您在图像的边框周围进行类似的检查，您会看到应该在那里的条目......没有。看一下第三行第五列。你会认为Ng = 1and Nr = 1，但我们在矩阵中没有看到。

再举一个例子。为什么是M(Ng,Nr) = 4, Ng = 2, Nr = 4？好吧，看看每个像素都有一个2。在这里我们可以成功地捕获8像素邻唯一有效的位置是row=2, col=4，row=3, col=3，row=3, col=4，row=4, col=3，和row=4, col=4。通过应用我们已经看到的相同算法，您将看到对于每个位置，Nr = 4. 因此，我们看到了Ng = 2, Nr = 4 四次的这种组合，这就是为什么将位置设置为 4。但是，在 中row=3, col=4，这实际上是Nr = 5，因为在该中心的附近有五个2。这就是为什么你看到Ng = 2, Nr = 5, M(Ng,Nr) = 1.

例如，让我们做一个位置。让我们2在矩阵 ( row=3, col=3)中间轻拍：

1. Ng = 2
2. 什么是有效的相邻像素？ 1, 1, 2, 0, 2, 3, 2, 2（省略中心）
3. 数一数等于 2 的像素数，一共有四个，所以 Nr = 4
4. M(Ng,Nr) += 1. 取Ng = 2，Nr = 4并将该点加 1。

如果您对具有 的其他有效位置执行此操作，则每次2都会看到，Nr = 4但第三行和第四列除外，其中Nr = 5。

那么我们如何在 MATLAB 中实现它呢？您可以做的是im2col将每个有效邻域转换为列。我还要做的是提取每个街区的中心。这实际上是矩阵的中间行。然后我们将计算出每个邻域有多少像素等于中心，将它们相加，这将确定我们的Nr值。这些Ng值将是中间行值本身。一旦我们这样做，我们就可以根据这些值计算直方图，就像算法如何获得我们的矩阵一样。换句话说，尝试这样做：

% // 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)

...这是我们的矩阵：

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

如果你检查这个，你会看到它与Q.

奖金

如果您希望能够在一般情况下适应算法，在您指定d和 的地方a，我们可以简单地遵循您的文本指南。对于每个邻域，您会发现中心像素与所有其他像素之间的差异。您计算<= a任何正整数的像素数d。请注意，这将创建一个2*d + 1 x 2*d + 1我们需要检查的邻域。我们也可以把它变成一个函数。无需再费周折：

%// 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

我们可以通过以下方式使用示例矩阵重新创建上面显示的场景：

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 是这样的：

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

希望这可以帮助！