两幅图像的互信息和联合熵 - MATLAB

Posted 2023-03-15

【中文标题】两幅图像的互信息和联合熵 - MATLAB

【英文标题】：Mutual information and joint entropy of two images - MATLAB

【发布时间】：2014-07-04 16:15:38

【问题描述】：

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

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

【问题讨论】：

如果您的图像强度在[0,1] 的范围内，那么您仍然可以使用我在下面编写的代码，但您需要确保强度是整数。如果您有 8 位图像，请在完成我所做的之前使用 im2uint8。

另外，如果您的图像是纯二进制的，则不要使用im2uint8 进行转换，因为这会浪费空间。您可以按原样使用代码。将其保留为具有 4 个联合箱的直方图，而不是 256 x 256。

非常感谢您的帮助。我最终将我的图像乘以 100，然后将数字四舍五入 + 1 -> Im1 = round(100*Im0+1) 这样我修复了所有错误 1）代码需要整数 2）它需要正值。 :D

如果您的图像是 8 位图像，我不建议您这样做，因为会丢失 bin 计数。我仍然会使用im2uint8 转换您的图像。这会将您的图像转换为[0,255]。此外，如果您查看 MATLAB 如何实现 entropy 命令，他们也会使用 im2uint8。您需要与 entropy 使用的 bin 数量以及我们讨论的联合熵保持一致。

【参考方案1】：

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

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

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

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

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

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

换句话说：

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

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

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

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

因此，联合熵可以计算为：

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

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

如果我们有浮点数据怎么办？

到目前为止，我们假设您处理的图像具有整数值的强度。如果我们有浮点数据怎么办？ 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));

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

一边

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

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;

---

希望这会有所帮助！

【讨论】：

使用 for 循环的不错选择。

这几乎是垃圾评论，但不得不说：这是我在 SO 上遇到的最佳答案之一。它简明扼要地解释了背景、代码背后的基本原理，并从简单（整数）到更难的实际情况（实际值）。此外，使用accumarray 的最简洁示例之一。

@ArunIRC 非常感谢你 :) 这实际上是我开始回答问题时的早期答案之一。没想到会这么受欢迎。实际上，您现在看到的当前版本与初始版本相比有很多修改。我只解决了整数情况，只解决了联合直方图。我随着时间的推移对其进行了编辑，并包括了浮点情况，最后计算了互信息。无论哪种方式，感谢您的评论:) 非常感谢。

@Keivan 感谢您的评论。实际上，有一个我没有发现的错误在您的示例中变得很明显。我刚刚修好了。现在运行它，它应该是正确的。基本上，我没有正确地创建两个信号之间的概率分布。但是，要回答您的问题，两个信号之间的熵为 0，这意味着它们之间的互信息也为 0。这是因为所有值仅映射到 1 个 bin，因此该值出现的概率始终为 1 ，因此其中的 log2 为 0。

@Keivan 您是否使用假设双精度的代码？您需要使用适合数据类型的代码。