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

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

Image 1 = X Image 2 = Y 

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

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

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

这里是我想要find联合熵的图像的一个例子:

 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 

要计算联合熵,需要计算两幅图像之间的联合直方图。 联合直方图基本上与正常的一维直方图相同,但是第一维度logging第一图像的强度和第二图像的第二维度logging强度。 这与通常所说的共生matrix非常相似。 在联合直方图中的位置(i,j) ,它告诉你我们遇到了多less强度值,其中第一个图像的强度为j ,第二个图像的强度为j

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

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

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

有人会倾向于用for循环来做这个,但是正如众所周知的那样,循环是非常缓慢的,如果可能的话应该避免。 但是,您可以通过以下方式轻松地在MATLAB中执行此操作, 而无需使用循环。 假设im1im2是您想要比较的第一张和第二张图片。 我们可以做的是将im1im2转换成vector。 然后我们可以使用accumarray来帮助我们计算联合直方图。 accumarray是MATLAB中最强大的function之一。 你可以把它想象成一个缩小的MapReduce范例。 简而言之,每个数据input都有一个键和一个相关的值。 accumarray的目标是将所有属于同一个键的值进行分类并对所有这些值进行一些操作。 在我们的例子中,“关键”将是强度值,而值本身是每个强度值的值1 。 然后,我们希望所有映射到同一个bin的值1 相加 ,这正是我们计算直方图的方式。 accumarray的默认行为是添加所有这些值。 具体而言, accumarray的输出将是一个数组,每个位置计算映射到该键的所有值的总和。 例如,第一个位置将是映射到1的键的所有值的总和,第二个位置将是映射到2的键的所有值的总和等等。

但是,对于联合直方图,您想要确定哪些值映射到(i,j)的相同强度对,因此这里的键将是一对二维坐标。 如此,在两个图像之间共享的具有在第一图像中的强度i和在第二图像中的相同空间位置中的 j任何强度都是相同的密钥。 因此,在2D情况下, accumarray的输出将是一个2Dmatrix,其中每个元素(i,j)包含映射到键(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 ,第一个input是键,第二个input是值。 accumarray是,如果每个键都有相同的值,你可以简单地给第二个input分配一个常量,这就是我所做的,它是1 。 通常,这是一个与第一个input具有相同行数的数组。 另外,请特别注意前两行。 在你的图像中不可避免地会有一个0的强度,但是因为MATLAB开始索引为1 ,所以我们需要将这两个数组偏移1

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

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

注意我search的是联合直方图,而不是联合概率matrix。 这是因为联合直方图由整数组成,而联合概率matrix将在01之间。 由于这种划分,我想避免将这个matrix中的任何条目与0进行比较,因为数值舍入和不稳定性。 以上也将把我们的联合概率matrix转换成一个一维向量的叠加向量,这很好。

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

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

如果我对计算MATLAB中图像的熵的理解是正确的,那么它应该计算256档上的直方图/概率分布,所以你可以肯定在这里用刚刚计算的联合熵来使用这个函数。

如果我们有浮点数据呢?

到目前为止,我们已经假定你所处理的图像具有整数值的强度。 如果我们有浮点数据呢? accumarray假设你正试图用整数索引到输出数组中,但是我们仍然可以通过这个小小的碰撞来完成我们想要的东西。 你会做的只是分配两个图像的每个浮点值有一个唯一的ID 。 因此,您将使用这些ID的accumarray代替。 为了便于这个ID分配,使用unique – 特别是函数的第三个输出。 你会把每一个图像,把它们放在unique ,使这些指标被input到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 ,我们直接为这些variables分配unique的第三个输出,然后使用之前计算的相同的联合熵代码。 我们也不必像以前那样将variables加1,因为unique从1开始分配ID。

在旁边

您可以使用联合概率matrix实际计算每个图像的直方图或概率分布。 如果要计算第一幅图像的直方图/概率分布,则只需对每行累积所有列即可。 为了做第二个图像,你可以简单地累积每一列的所有行。 因此,你可以这样做:

 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) / numel(histogramImage1); %// Compute the entropy entropy1 = -sum(prob1NoZero.*log2(prob1NoZero)); %// Repeat for the second image indNoZero = histogramImage2 ~= 0; prob2NoZero = histogramImage2(indNoZero) / numel(histogramImage2); entropy2 = -sum(prob2NoZero.*log2(prob2NoZero)); %// Now compute mutual information mutualInformation = entropy1 + entropy2 - jointEntropy; 

希望这可以帮助!