29

我有两张黑白图像,我需要计算互信息。

Image 1 = X 
Image 2 = Y

我知道互信息可以定义为:

MI = entropy(X) + entropy(Y) - JointEntropy(X,Y)

MATLAB 已经内置了计算熵但不计算联合熵的函数。我想真正的问题是:如何计算两张图像的联合熵?

这是我想找到联合熵的图像示例:

X =

0    0    0    0    0    0
0    0    1    1    0    0
0    0    1    1    0    0
0    0    0    0    0    0
0    0    0    0    0    0

Y =

0    0    0    0    0    0 
0    0    0.38 0.82 0.38 0.04 
0    0    0.32 0.82 0.68 0.17
0    0    0.04 0.14 0.11 0 
0    0    0    0    0    0
4

1 回答 1

67

要计算联合熵,需要计算两幅图像之间的联合直方图。联合直方图与正常的一维直方图基本相同,但第一维记录第一张图像的强度,第二维记录第二张图像的强度。这与通常所说的共现矩阵非常相似。在联合直方图中的位置(i,j),它告诉您我们遇到了多少强度值,i在第一张图像中具有强度,j在第二张图像中具有强度。

重要的是,它记录了我们在相同对应位置看到这对强度的次数。例如,如果我们有一个联合直方图计数(7,3) = 2,这意味着当我们扫描两幅图像时,当我们遇到强度 时7,在第二张图像的相同对应位置,我们3总共遇到强度2

构建联合直方图非常简单。

  1. 首先,创建一个256 x 256矩阵(假设您的图像是无符号 8 位整数)并将它们初始化为全零。此外,您需要确保两个图像的大小(宽度和高度)相同。
  2. 完成此操作后,请查看每个图像的第一个像素,我们将其表示为左上角。具体来说,查看该位置的第一张和第二张图像的强度。第一张图像的强度将作为行,而第二张图像的强度将作为列。
  3. 在矩阵中找到这个位置,并将矩阵中的这个位置增加1.
  4. 对图像中的其余位置重复此操作。
  5. 完成后,将所有条目除以任一图像中的元素总数(记住它们的大小应该相同)。这将为我们提供两幅图像之间的联合概率分布。

人们倾向于使用for循环来执行此操作,但众所周知,for循环速度非常慢,应尽可能避免。但是,您可以通过以下方式在 MATLAB 中轻松执行此操作,而无需循环。假设im1im2是您要比较的第一张和第二张图像。我们能做的就是将im1和转换im2成向量。然后我们可以用它accumarray来帮助我们计算联合直方图。 accumarray是 MATLAB 中最强大的函数之一。您可以将其视为一个微型 MapReduce 范例。简单地说,每个数据输入都有一个键和一个关联的值。的目标accumarray就是对属于同一个key的所有值进行bin,并对所有这些值做一些操作。在我们的例子中,“关键”是强度值,而这些值本身就是1每个强度值的值。然后,我们希望将该映射的所有值加到1同一个 bin 中,这正是我们计算直方图的方式。的默认行为accumarray是添加所有这些值。具体来说, 的输出accumarray将是一个数组,其中每个位置计算映射到该键的所有值的总和。例如,第一个位置是映射到键 1 的所有值的总和,第二个位置是映射到键 2 的所有值的总和,依此类推。

但是,对于联合直方图,您想弄清楚哪些值映射到相同的强度对(i,j),因此这里的键将是一对 2D 坐标。这样,i在第一图像和j第二图像中在两个图像之间共享的相同空间位置中具有强度的任何强度都转到相同的键。因此,在 2D 情况下, 的输出accumarray将是一个 2D 矩阵,其中每个元素(i,j)包含映射到 key 的所有值的总和(i,j),类似于前面提到的 1D 情况,这正是我们所追求的。

换句话说:

indrow = double(im1(:)) + 1;
indcol = double(im2(:)) + 1; %// Should be the same size as indrow
jointHistogram = accumarray([indrow indcol], 1);
jointProb = jointHistogram / numel(indrow);

使用accumarray时,第一个输入是键,第二个输入是值。需要注意的accumarray是,如果每个键具有相同的值,您可以简单地为第二个输入分配一个常量,这就是我所做的,它是1. 通常,这是一个与第一个输入具有相同行数的数组。另外,请特别注意前两行。您的图像中不可避免地会有一个强度0,但由于 MATLAB 从 开始索引1,我们需要将两个数组偏移1

现在我们有了联合直方图,计算联合熵真的很简单。它类似于 1D 中的熵,只是现在我们只是对整个联合概率矩阵求和。请记住,您的联合直方图很可能有很多0条目。我们需要确保跳过这些,否则log2操作将是未定义的。现在让我们摆脱任何零条目:

indNoZero = jointHistogram ~= 0;
jointProb1DNoZero = jointProb(indNoZero);

请注意,我搜索的是联合直方图而不是联合概率矩阵。这是因为联合直方图由整数组成,而联合概率矩阵位于0和之间1。由于除法,我想避免将此矩阵中的任何条目与0数字舍入和不稳定性进行比较。以上还将我们的联合概率矩阵转换为堆叠的一维向量,这很好。

因此,联合熵可以计算为:

jointEntropy = -sum(jointProb1DNoZero.*log2(jointProb1DNoZero));

如果我对在 MATLAB 中计算图像熵的理解是正确的,它应该计算 bin 上的直方图/概率分布256,因此您当然可以在此处使用该函数和刚刚计算的联合熵。

如果我们有浮点数据呢?

到目前为止,我们假设您处理的图像具有整数值的强度。如果我们有浮点数据怎么办? accumarray假设您正在尝试使用整数对输出数组进行索引,但我们仍然可以通过这个小小的颠簸来完成我们想要的。您要做的只是将两个图像中的每个浮点值分配为具有唯一的 ID。因此,您将accumarray改为使用这些 ID。为了便于分配此 ID,请使用unique- 特别是函数的第三个输出。您将拍摄每张图像,将它们放入unique并将这些索引输入到accumarray. 换句话说,改为这样做:

[~,~,indrow] = unique(im1(:)); %// Change here
[~,~,indcol] = unique(im2(:)); %// Change here

%// Same code
jointHistogram = accumarray([indrow indcol], 1);
jointProb = jointHistogram / numel(indrow);
indNoZero = jointHistogram ~= 0;
jointProb1DNoZero = jointProb(indNoZero);
jointEntropy = -sum(jointProb1DNoZero.*log2(jointProb1DNoZero));

请注意,使用indrowindcol,我们直接将 的第三个输出分配unique给这些变量,然后使用我们之前计算的相同联合熵代码。我们也不必像以前那样将变量偏移 1,因为将从1 开始unique分配 ID 。

在旁边

您实际上可以使用联合概率矩阵单独计算每个图像的直方图或概率分布。如果您想计算第一张图像的直方图/概率分布,您只需累积每行的所有列。要为第二个图像执行此操作,您只需累积每列的所有行。因此,您可以这样做:

histogramImage1 = sum(jointHistogram, 1);
histogramImage2 = sum(jointHistogram, 2);

之后,您可以自己计算这两者的熵。要仔细检查,请确保将两者都转换为 PDF,然后使用标准方程计算熵(如上)。


我如何最终计算互信息?

要最终计算互信息,您将需要两个图像的熵。您可以使用 MATLAB 的内置entropy函数,但这假设有 256 个唯一级别。您可能希望将其应用于存在N不同级别而不是 256 的情况,因此您可以使用我们上面对联合直方图所做的操作,然后在上面的备用代码中计算每个图像的直方图,然后计算熵对于每个图像。您只需重复联合使用的熵计算,但将其单独应用于每个图像:

%// Find non-zero elements for first image's histogram
indNoZero = histogramImage1 ~= 0;

%// Extract them out and get the probabilities
prob1NoZero = histogramImage1(indNoZero);
prob1NoZero = prob1NoZero / sum(prob1NoZero);

%// Compute the entropy
entropy1 = -sum(prob1NoZero.*log2(prob1NoZero));

%// Repeat for the second image
indNoZero = histogramImage2 ~= 0;
prob2NoZero = histogramImage2(indNoZero);
prob2NoZero = prob2NoZero / sum(prob2NoZero);
entropy2 = -sum(prob2NoZero.*log2(prob2NoZero));

%// Now compute mutual information
mutualInformation = entropy1 + entropy2 - jointEntropy;

希望这可以帮助!

于 2014-05-16T03:05:36.690 回答