如何从一个清晰的例子中计算二维图像中的吉布斯能量 [英] How to compute Gibbs energy in 2D image from a clear example
问题描述
我有一个关于矩阵的有趣问题。在吉布斯分布中，吉布斯能量 U（x）可以计算为
I have an interesting question about matrix. In Gibbs distribution, a Gibbs energy U(x) can be compute as
这是所有可能派系C（右图）上的集团潜力Vc（x）的总和。集团c被定义为S中的站点子集（x-蓝色像素的邻域是左图中的黄色像素的邻居），其中每对不同的站点是邻居，除了单站点集团
which is a sum of clique potentials Vc(x) over all possible cliques C (right image). A clique c is defined as a subset of sites in S (are neighborhood of x- blue color pixels are neighbor of yellow pixel in left figure) in which every pair of distinct sites are neighbors, except for single-site cliques
其中 V（x）计算如下：
< img src =https://i.stack.imgur.com/8ojSX.jpgalt =在此处输入图像说明>
我的目标是如何计算 U（x）。在上面的示例中， x = {1,2} 。但是，我有一些卡在图像角落的像素，只有少于8个邻居（正常情况下像素通常有8个相邻像素）。为了解决这个问题，我通过带有掩码的卷积图像在边角处添加了零。但我认为它可能会对U（x）产生影响。从上面的定义，当我考虑图像角落处的像素（灰色）时，你能帮我找到U（x = 1）和U（x = 2）吗？我尝试通过以下解决方案解决它，但我认为我的解决方案太长了，我不确定是否正确。
My goal is how to compute U(x). In my above example, x={1,2}. However, I have some stuck at pixel at corner of image that only has less than 8 neighbors (a pixel often has 8 neighbor pixels in normal case). To solve it, I added zeros at corner of edge by convolution image with an mask. But I think it maybe affect to U(x). From above definition, could you help me find U(x=1) and U(x=2) when i consider pixels at corner of image (gray color)? I tried to solve it by below solution but I think my solution is so long and I am not sure is whether correct or not.
这就是我做的事情
Imlabel =[ 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 2 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 2 2 2 1 1 1 1 1;
1 1 2 2 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 2 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1];
mask=[ 0 0 0;0 1 0;0 0 0];
Imlabel=conv2(Imlabel,mask); % To solve pixel in corner
num_class=2; % x={1,2}
beta=1;
for label=1:num_class
for i=2:size(Imlabel,1)-1
for j=2:size(Imlabel,2)-1
sum_V=0;
%North, south, east and west neighbors
if(label==Imlabel(i-1,j)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i,j+1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i+1,j)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i,j-1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
%% Diagonal pixels
if(label==Imlabel(i-1,j-1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i-1,j+1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i+1,j+1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i+1,j-1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
%% Save
U(i-1,j-1,label)=sum_V; %% Because Index is extended before
end
end
end
更新：标准化
目前，基于我对此处。我计算了归一化项以及gibbs分布为
Update: Normalization Currently, based on my understand at here. I computed normalized term as well as gibbs distribution as
P=zeros(size(U));
Z=sum(exp(-U),3);
for i=1:num_class
P(:,:,i)=exp(-U(:,:,i))./Z;
end
推荐答案
你的解释也是有道理的作为代码。但是，你可以通过聪明地想要使用的各种函数来完成这个矢量化。
Your explanation makes sense as well as the code. However, you can totally get this vectorized by being smart about the kinds of functions you want to use.
对于感兴趣的读者，给出你的源代码，你想要什么do是以下（伪代码）：
For the interested readers, given your source code, what you want to do is the following (in pseudocode):
-
创建矩阵
U那是尺寸（Imlabel，1）x尺寸（Imlabel，2）x num_class大。该矩阵的每个切片将确定该切片中每个坐标（i，j）处的8个相邻像素的吉布斯能量成本。
Create a matrix
Uthat issize(Imlabel,1) x size(Imlabel,2) x num_classlarge. Each slice of this matrix will determine which the Gibbs energy costs for the 8 neighbouring pixels at each coordinate(i,j)in this slice.
对于每个类别标签 x ...
For each class label x...
a。对于图像中的每个像素（i，j） ...
a. For each pixel in the image (i,j)...
-
找到在
（i，j）的8像素邻域中定义的等于x的位置将它们设置为beta
Find those locations defined in an 8-pixel neighbourhood of
(i,j)that are equal toxand set them tobeta
查找在<$ c $的8像素邻域中定义的位置c>（i，j）不等于 x 并将它们设置为 -beta
Find those locations defined in an 8-pixel neighbourhood of (i,j) that are not equal to x and set them to -beta
计算所有这些 beta 和 -beta的总和此8像素区域内的值，并将其设置为 U（i，j，x）
Compute the sum of all of these beta and -beta values within this 8-pixel neighbourhood and set this at location U(i,j,x)
因此，我要做的是创建一个 3D矩阵的成本。 ..让我们称之为 C ，其中每个切片 y 的尺寸与尺寸相同（ Imlabel，1）x size（Imlabel，2）但是切片中的每个位置（i，j）将是 beta if Imlabel（i，j）== y 和 -beta 否则。执行此操作后，您可以对此矩阵执行ND卷积，但使用2D内核，以便您可以分别计算每个切片的8像素邻域的总和。
Therefore, what I would do is create a 3D matrix of costs... let's call this C, where each slice y has the same dimensions as size(Imlabel,1) x size(Imlabel,2) but each location (i,j) in a slice will be beta if Imlabel(i,j) == y and -beta otherwise. Once you do this, you can perform an N-D convolution of this matrix, but using a 2D kernel so that you can compute the summation of the 8-pixel neighbourhoods of each slice separately.
今天的魔术菜单包括 bsxfun ， convn 以及 置换 以增加底池效果。
Today's magic menu consists of bsxfun, convn and a side-order of permute to sweeten the pot.
要获得成本矩阵 C ，您可以执行此操作。这假设 Imlabel 沿着边界没有用零填充：
To get the cost matrix C, you can do this. This is assuming Imlabel is not padded with zeros along the boundaries:
Imlabel =[ 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 2 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 2 2 2 1 1 1 1 1;
1 1 2 2 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 2 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1];
C = double(bsxfun(@eq, Imlabel, permute(1:num_class, [1 3 2])));
C(C == 0) = -beta;
C(C == 1) = beta;
permute 函数用于创建单个3D矢量，从1到最多为您定义的类。之所以需要这样做是因为 bsxfun 执行所谓的广播。在这种情况下会发生什么，你的矩阵 Imlabel ，使用3D矢量将与 eq 或等于功能。此成本矩阵的每个切片都会为您提供与相关标签相等的位置，或者不是。
The permute function is used to create a single 3D vector that goes from 1 up to as many classes as you have defined. The reason why this is required is because bsxfun does what is known as broadcasting. What will happen in this case is given your matrix Imlabel, using the 3D vector will create a 3D cost matrix in conjunction with the eq or equals function. Each slice of this cost matrix will give you the locations of where it is either equal to the label in question, or it isn't.
您可以验证我们的内容从 bsxfun 输出的位置矩阵是正确的：
You can verify that what we have for the locations matrix that is output from bsxfun is correct:
>> C = double(bsxfun(@eq, Imlabel, permute(1:num_class, [1 3 2])))
C(:,:,1) =
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 0 1 1
1 1 1 1 1 1 1 1 1 1
1 1 0 0 0 1 1 1 1 1
1 1 0 0 1 1 1 1 1 1
1 1 1 1 1 1 1 1 0 1
1 1 1 1 1 1 1 1 1 1
1 1 1 1 1 1 1 1 1 1
C(:,:,2) =
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 1 0 0
0 0 0 0 0 0 0 0 0 0
0 0 1 1 1 0 0 0 0 0
0 0 1 1 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0
0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0
如您所见，在每个切片中，我们可以看到哪些位置共享该切片枚举的标签以及那些不共享相同标签的标签。一旦我们有这些位置，我们需要将其转换为 double ，以便我们可以将其修改为 -beta 是一个位置，它不等于特定切片的标签， beta 是。这是通过 bsxfun 调用之后的两个赋值语句完成的。
As you can see, for each slice, we can see which locations shared that label that is enumerated by that slice and those that don't share the same label. Once we have these locations, we need to cast this to double so that we can modify this to be a cost matrix where -beta is a location that isn't equal to a label for a particular slice and beta is. That is done by the two assignment statements after the bsxfun call.
一旦你有了这个成本矩阵，就可以像以前一样填充边界，但要确保边界设置为的成本-beta 就像你在代码中所做的那样：
Once you have this cost matrix, you can pad the boundaries like you did before, but make sure that the boundaries are set to a cost of -beta as you have done in your code:
mask = zeros(3,3); mask(2,2) = 1;
Cpad = convn(C, mask);
Cpad(Cpad == 0) = -beta;
第一行代码表示 [0 0 0; 0 1 0; 0 0 0] 您在代码中定义的掩码。现在要单独填充每个切片以使其具有零边框，我们可以使用 convn 来帮助我们这样做。
The first line of code denotes the [0 0 0; 0 1 0; 0 0 0] mask that you defined in your code. Now to pad each slice independently to have a border of zeroes, we can use convn to help us do that.
现在我们有了这个设置，你所要做的就是每片独立地计算Gibbs的能量：
Now that we have that set up, all you have to do is perform the calculation of Gibbs' energy independently per slice:
mask2 = ones(3,3); mask2(2,2) = 0;
out = convn(Cpad, mask2, 'valid');
第一行代码表示一个掩码，除了中间值，每个值为1，即0它是一个3 x 3的面具。你想要这样做的原因是因为它取代了你的双中的循环逻辑。你正在做的基本上是一个卷积，你将所有8个邻居加在一起，除了中间本身。这可以通过使用除中心像素之外的所有1的掩码来实现，并将其设置为0. if / else 语句，用于
The first line of code denotes a mask where every value is 1 except for the middle, which is 0 and it is a 3 x 3 mask. The reason why you want to do this is because this replaces the if/else statements in your double for loop logic. What you are doing is essentially a convolution where you are adding up all of the 8 neighbours except for the middle itself. This can be accomplished by using a mask that has all 1s except for the central pixel and you set that to 0.
接下来，我们使用 convn 独立计算每个切片的吉布斯能量，但使用'有效'标志，以便我们可以删除我们包含的零填充边框开头。
Next, we use convn to compute the Gibbs' energy per slice independently, but using the 'valid' flag so that we can remove the zero-padded border we included at the beginning.
现在，如果你看一下 out ，它会与你用 U ，但由于您将 U 编入索引的方式略有改变。不过，上面输出 out ，你可以验证我们所拥有的是正确的，这将很好地处理边界情况：
Now, if you take a look at out, it compares with what you have with U, but there is a slight shift because of the way you are indexing into U. Nevertheless, with the above output of out, you can verify that what we have is correct and this will handle the border cases nicely:
>> out
out(:,:,1) =
-2 2 2 2 2 2 2 2 2 -2
2 8 8 8 8 8 6 6 6 2
2 8 8 8 8 8 6 8 6 2
2 6 4 2 4 6 6 6 6 2
2 4 2 0 4 6 8 8 8 2
2 4 2 0 2 6 8 6 6 0
2 6 4 4 6 8 8 6 8 0
2 8 8 8 8 8 8 6 6 0
-2 2 2 2 2 2 2 2 2 -2
out(:,:,2) =
-8 -8 -8 -8 -8 -8 -8 -8 -8 -8
-8 -8 -8 -8 -8 -8 -6 -6 -6 -8
-8 -8 -8 -8 -8 -8 -6 -8 -6 -8
-8 -6 -4 -2 -4 -6 -6 -6 -6 -8
-8 -4 -2 0 -4 -6 -8 -8 -8 -8
-8 -4 -2 0 -2 -6 -8 -6 -6 -6
-8 -6 -4 -4 -6 -8 -8 -6 -8 -6
-8 -8 -8 -8 -8 -8 -8 -6 -6 -6
-8 -8 -8 -8 -8 -8 -8 -8 -8 -8
例如，如果我们看一下第一片的左上角，我们看到东，南，东南角都有相同的标签1，因此 beta = 1 ，我们的总和将是3，但是我们没有考虑其他5个方向，并在你的代码，你将它们设置为 -beta ，所以 3 - 5 = -2 。
For example, if we took a look at the top left corner of the first slice, we see that the east, south and south east corners have the same label of 1, and so since beta = 1, our sum would be 3, but then there are the other 5 directions that we didn't consider, and in your code, you set those to -beta and so 3 - 5 = -2.
让我们看看第6行第4列，让我们来看看第二个标签或切片。这意味着任何等于标签2的基本方向，我们应该将这些与 beta / 1 相加，而不等于标签2的任何东西都是 -beta / -1 。正如你所看到的，正好有4个1和4个标签的标签为2，所以这些应该取消并给我们一个0。
Let's also take a look at the 6th row, 4th column and let's take a look at the second label or slice. This means that any cardinal directions that are equal to label 2, we should sum these with beta / 1 while anything that isn't equal to label 2 is -beta / -1. As you can see, there are exactly 4 labels of 1 and 4 labels of 2, and so these should cancel out and give us a 0.
你可以验证那是什么我们为此矩阵中的所有其他位置做了正确的计算。
You can verify that what we did for all of the other locations in this matrix result in the right calculations.
完整代码只是：
Imlabel =[ 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 2 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 2 2 2 1 1 1 1 1;
1 1 2 2 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 2 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1];
C = double(bsxfun(@eq, Imlabel, permute(1:num_class, [1 3 2])));
C(C == 0) = -beta;
C(C == 1) = beta;
mask = zeros(3,3); mask(2,2) = 1;
Cpad = convn(C, mask);
Cpad(Cpad == 0) = -beta;
mask2 = ones(3,3); mask2(2,2) = 0;
out = convn(Cpad, mask2, 'valid');
在时间的情况下，我们可以检查看看上述方法与原始循环代码的比较速度有多快。我使用 timeit 帮助促进这一点。
In the case of timing, we can check to see how fast the above approach is compared to your original loop code. I used timeit to help facilitate that.
这是我设置的脚本，用于设置所有相关变量，并测量代码和我的代码所需的时间：
Here's the script I set up that sets up all of the relevant variables, and measures the time taken with your code and mine:
function test_clique_timing
% Setup
Imlabel_orig =[ 1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 2 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1;
1 1 2 2 2 1 1 1 1 1;
1 1 2 2 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 2 1;
1 1 1 1 1 1 1 1 1 1;
1 1 1 1 1 1 1 1 1 1];
num_class=2; % x={1,2}
beta = 1;
function test_orig
mask=[ 0 0 0;0 1 0;0 0 0];
Imlabel=conv2(Imlabel_orig,mask); % To solve pixel in corner
beta=1;
for label=1:num_class
for i=2:size(Imlabel,1)-1
for j=2:size(Imlabel,2)-1
sum_V=0;
%North, south, east and west neighbors
if(label==Imlabel(i-1,j)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i,j+1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i+1,j)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i,j-1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
%% Diagonal pixels
if(label==Imlabel(i-1,j-1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i-1,j+1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i+1,j+1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
if(label==Imlabel(i+1,j-1)) sum_V=sum_V+beta;
else sum_V=sum_V-beta;end
%% Save
U(i-1,j-1,label)=sum_V; %% Because Index is extended before
end
end
end
end
function test_conv
C = double(bsxfun(@eq, Imlabel_orig, permute(1:num_class, [1 3 2])));
C(C == 0) = -beta;
C(C == 1) = beta;
mask = zeros(3,3); mask(2,2) = 1;
Cpad = convn(C, mask);
Cpad(Cpad == 0) = -beta;
mask2 = ones(3,3); mask2(2,2) = 0;
out = convn(Cpad, mask2, 'valid');
end
time1 = timeit(@test_orig);
time2 = timeit(@test_conv);
fprintf('Loop code time: %f seconds\n', time1);
fprintf('Vectorized code time: %f seconds\n', time2);
fprintf('Speedup factor: %f', time1/time2);
end
运行上面的脚本，我有时会得到这些：
Running the above script, I get these for times:
Loop code time: 0.000049 seconds
Vectorized code time: 0.000060 seconds
Speedup factor: 0.816667
这是在几秒钟内完成的，我是在Mac OS Mavericks 10.10.5下使用MATLAB R2013a完成的。加速看起来并不那么好。事实上它是一个小于1的因素，这很糟糕。但是，您所展示的是如此小的数据集。我们应该尝试使用更大的数据集，看看它是否仍然存在。让我们制作矩阵2500 x 2500，有10个标签。我将用以下代码替换 Imlabel 矩阵：
This is done in seconds and I did this using MATLAB R2013a under Mac OS Mavericks 10.10.5. The speedup doesn't look that great. In fact it's a factor less than 1, which is terrible. However, what you have shown is such a small dataset. We should try with a larger dataset and see if this still holds. Let's make the matrix 2500 x 2500, with 10 class labels. I'm going to replace the Imlabel matrix with:
rng(123); %// Set seed for reproducibility
num_class = 10;
Imlabel_orig = randi(10,2500,2500);
这会在2500 x 2500矩阵中生成1到10的随机整数。这样做并再次运行代码给出：
This produces random integers from 1 to 10 in a 2500 x 2500 matrix. Doing this and running the code again gives:
Loop code time: 17.553669 seconds
Vectorized code time: 3.321950 seconds
Speedup factor: 5.284146
是的......那好多了。在2500 x 2500矩阵的10个标签上，计算需要大约17.5秒，而带有 bsxfun / convn / permute 的矢量化代码优于原始循环代码的因子差不多是5.2倍。
Yeah... that's much better. On the 2500 x 2500 matrix of 10 labels, it took roughly 17.5 seconds to compute, whereas the the vectorized code with bsxfun / convn / permute outperforms your original loop code by a factor of almost 5.2x.
至于影响你的吉布斯能量计算的归一化项，你现在有：
As for the normalization term to factor into your Gibbs' energy calculation, you currently have this:
P=zeros(size(U));
Z=sum(exp(-U),3);
for i=1:num_class
P(:,:,i)=exp(-U(:,:,i))./Z;
end
如果你利用 bsxfun ，你可以这样做：
If you take advantage of bsxfun, you can simply do this:
Z = sum(exp(-out),3)); %// out is U in my code
P = bsxfun(@rdivide, exp(-out), Z);
这实现了代码所做的一样。对于 P 中的每个切片，我们找到 exp 并将切片作为输入否定，并将每个项除以ž。对于 bsxfun ， Z 矩阵的复制次数与 out 并且元素分区与循环代码非常相似。
This achieves the same thing your code does. For each slice in P, we find the exp and negate the slice as input, and divide each term by Z. For bsxfun, the Z matrix gets replicated for as many slices as there are in out and does the element-wise division much like the loop code does.
这篇关于如何从一个清晰的例子中计算二维图像中的吉布斯能量的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
![< img src =https://i.stack.imgur.com/8ojSX.jpgalt =在此处输入图像说明>](\"https://i.stack.imgur.com/8ojSX.jpg\"){kind=link}