两个图像的互信息和联合熵 - MATLAB [英] Mutual information and joint entropy of two images - MATLAB

问题描述

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

I have two black and white images and I need to calculate the mutual information.

Image 1 = X 
Image 2 = Y

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

I know that the mutual information can be defined as:

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

MATLAB已经有内置函数来计算熵而不是计算联合熵。我想真正的问题是：我如何计算两幅图像的联合熵？

MATLAB already has built-in functions to calculate the entropy but not to calculate the joint entropy. I guess the true question is: How do I calculate the joint entropy of two images?

以下是我想找到联合熵的图像示例：

Here is an example of the images I'd like to find the joint entropy of:

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

推荐答案

要计算联合熵，需要计算两幅图像之间的联合直方图。联合直方图与普通1D直方图基本相同，但第一维记录第一图像的强度，第二维记录第二图像的强度。这与通常所说的共现矩阵非常相似。在联合直方图中的位置（i，j），它会告诉您我们遇到的强度值 i 在第一张图片中，强度 j 在第二张图片中。

To calculate the joint entropy, you need to calculate the joint histogram between two images. The joint histogram is essentially the same as a normal 1D histogram but the first dimension logs intensities for the first image and the second dimension logs intensities for the second image. This is very similar to what is commonly referred to as a co-occurrence matrix. At location (i,j) in the joint histogram, it tells you how many intensity values we have encountered that have intensity i in the first image and intensity j in the second image.

重要的是，这记录了我们在相同的相应位置看到这对强度的次数。例如，如果我们的联合直方图计数（7,3）= 2 ，这意味着当我们扫描两个图像时，当我们遇到强度<$时c $ c> 7 ，在第二张图片的同一位置，我们遇到 3 的强度，总计 2 次。

What is important is that this logs how many times we have seen this pair of intensities at the same corresponding locations. For example, if we have a joint histogram count of (7,3) = 2, this means that when we were scanning both images, when we encountered the intensity of 7, at the same corresponding location in the second image, we encountered the intensity of 3 for a total of 2 times.

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

Constructing a joint histogram is very simple to do.

1. 首先，创建一个 256 x 256 矩阵（假设您的图像是无符号的8 -bit integer）并将它们初始化为全零。此外，您需要确保两个图像的大小（宽度和高度）相同。


2. 一旦你这样做，看看每个图像的第一个像素，我们将其表示为左上角。具体来说，请查看此位置的第一张和第二张图像的强度。第一图像的强度将用作行，而第二图像的强度将用作列。


3. 在矩阵中找到此位置，并通过 1 在矩阵中增加此点。


4. 对图像中的其余位置重复此操作。


5. 完成后，将所有条目除以中的元素总数任何一个图像（记住它们应该是相同的大小）。这将为我们提供两个图像之间的联合概率分布。

1. First, create a 256 x 256 matrix (assuming your image is unsigned 8-bit integer) and initialize them to all zeroes. Also, you need to make sure that both of your images are the same size (width and height).
2. Once you do that, take a look at the first pixel of each image, which we will denote as the top left corner. Specifically, take a look at the intensities for the first and second image at this location. The intensity of the first image will serve as the row while the intensity of the second image will serve as the column.
3. Find this location in the matrix and increment this spot in the matrix by 1.
4. Repeat this for the rest of the locations in your image.
5. After you're done, divide all entries by the total number of elements in either image (remember they should be the same size). This will give us the joint probability distribution between both images.

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

One would be inclined to do this with for loops, but as it is commonly known, for loops are notoriously slow and should be avoided if at all possible. However, you can easily do this in MATLAB in the following way without loops. Let's assume that im1 and im2 are the first and second images you want to compare to. What we can do is convert im1 and im2 into vectors. We can then use accumarray to help us compute the joint histogram. accumarray is one of the most powerful functions in MATLAB. You can think of it as a miniature MapReduce paradigm. Simply put, each data input has a key and an associated value. The goal of accumarray is to bin all of the values that belong to the same key and do some operation on all of these values. In our case, the "key" would be the intensity values, and the values themselves are the value of 1 for every intensity value. We would then want to add up all of the values of 1 that map to the same bin, which is exactly how we'd compute a histogram. The default behaviour for accumarray is to add all of these values. Specifically, the output of accumarray would be an array where each position computes the sum of all values that mapped to that key. For example, the first position would be the summation of all values that mapped to the key of 1, the second position would be the summation of all values that mapped to the key of 2 and so on.

然而，对于联合直方图，您想要找出哪些值映射到（i，j），所以这里的键是一对2D坐标。因此，第一张图片中强度为 i 的强度和第二张图片中 j 的强度在两个图像之间共享的相同空间位置转到相同的键。因此，在2D情况下， accumarray 的输出将是一个2D矩阵，其中每个元素（i，j）包含映射到键（i，j）的所有值的总和，类似于之前提到的1D情况，这正是我们所追求的。

However, for the joint histogram, you want to figure out which values map to the same intensity pair of (i,j), and so the keys here would be a pair of 2D coordinates. As such, any intensities that have an intensity of i in the first image and j in the second image in the same spatial location shared between the two images go to the same key. Therefore in the 2D case, the output of accumarray would be a 2D matrix where each element (i,j) contains the summation of all values that mapped to key (i,j), similar to the 1D case that was mentioned previously which is exactly what we are after.

换句话说：

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 。

With accumarray, the first input are the keys and the second input are the values. A note with accumarray is that if each key has the same value, you can simply assign a constant to the second input, which is what I've done and it's 1. In general, this is an array with the same number of rows as the first input. Also, take special note of the first two lines. There will inevitably be an intensity of 0 in your image, but because MATLAB starts indexing at 1, we need to offset both arrays by 1.

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

Now that we have the joint histogram, it's really simple to calculate the joint entropy. It is similar to the entropy in 1D, except now we are just summing over the entire joint probability matrix. Bear in mind that it will be very likely that your joint histogram will have many 0 entries. We need to make sure that we skip those or the log2 operation will be undefined. Let's get rid of any zero entries now:

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

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

Take notice that I searched the joint histogram instead of the joint probability matrix. This is because the joint histogram consists of whole numbers while the joint probability matrix will lie between 0 and 1. Because of the division, I want to avoid comparing any entries in this matrix with 0 due to numerical roundoff and instability. The above will also convert our joint probability matrix into a stacked 1D vector, which is fine.

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

As such, the joint entropy can be calculated as:

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

如果我对在MATLAB中计算图像的熵的理解是正确的，它应该计算直方图/概率分布在 256 bins，所以你当然可以在这里使用刚刚计算出来的联合熵。

If my understanding of calculating entropy for an image in MATLAB is correct, it should calculate the histogram / probability distribution over 256 bins, so you can certainly use that function here with the joint entropy that was just calculated.

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

So far, we have assumed that the images that you have dealt with have intensities that are integer-valued. What if we have floating point data? accumarray assumes that you are trying to index into the output array using integers, but we can still certainly accomplish what we want with this small bump in the road. What you would do is simply assign each floating point value in both images to have a unique ID. You would thus use accumarray with these IDs instead. To facilitate this ID assigning, use unique - specifically the third output from the function. You would take each of the images, put them into unique and make these the indices to be input into accumarray. In other words, do this instead:

[~,~,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 将从1 开始分配ID 。

Note that with indrow and indcol, we are directly assigning the third output of unique to these variables and then using the same joint entropy code that we computed earlier. We also don't have to offset the variables by 1 as we did previously because unique will assign IDs starting at 1.

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

You can actually calculate the histograms or probability distributions for each image individually using the joint probability matrix. If you wanted to calculate the histograms / probability distributions for the first image, you would simply accumulate all of the columns for each row. To do it for the second image, you would simply accumulate all of the rows for each column. As such, you can do:

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

之后，您可以自己计算这两者的熵。要仔细检查，请确保将这两个转换为PDF，然后使用标准公式计算熵（如上所述）。

After, you can calculate the entropy of both of these by yourself. To double check, make sure you turn both of these into PDFs, then compute the entropy using the standard equation (like above).

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

To finally compute Mutual Information, you're going to need the entropy of the two images. You can use MATLAB's built-in entropy function, but this assumes that there are 256 unique levels. You probably want to apply this for the case of there being N distinct levels instead of 256, and so you can use what we did above with the joint histogram, then computing the histograms for each image in the aside code above, and then computing the entropy for each image. You would simply repeat the entropy calculation that was used jointly, but apply it to each image individually:

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

---

希望这会有所帮助！

---

Hope this helps!

这篇关于两个图像的互信息和联合熵 - MATLAB的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！