matlab - 如何从一个清晰的例子中计算2D图像中的吉布斯能量

Question

我对矩阵有一个有趣的问题。在吉布斯分布中，吉布斯能量U（x）可计算为



这是所有可能的团C上团势Vc（x）的总和（右图）。集团c定义为S中站点的子集（左图中x-蓝色像素的邻居是黄色像素的邻居），其中每对不同的站点都是邻居，除了单站点集团



其中V（x）的计算如下：



我的目标是如何计算U（x）。在我上面的示例中，x = {1,2}。但是，我有一些卡在图像角上的像素中，该像素只有少于8个邻居（正常情况下一个像素通常有8个邻居像素）。为了解决这个问题，我在带遮罩的卷积图像的边缘角处添加了零。但我认为这可能会影响U（x）。根据上面的定义，当我考虑图像角点（灰色）时的像素时，可以帮我找到U（x = 1）和U（x = 2）吗？我试图通过下面的解决方案来解决它，但是我认为我的解决方案太长了，我不确定是否正确。

这就是我所做的

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

更新：规范化
目前，根据我对 here的了解。我计算归一化的术语以及吉布斯分布为

P=zeros(size(U));
Z=sum(exp(-U),3);
for i=1:num_class
   P(:,:,i)=exp(-U(:,:,i))./Z;
end

Answer 1

您的解释和代码一样有意义。但是，通过精通要使用的功能种类，可以完全实现矢量化。

对于感兴趣的读者，给定您的源代码，您想要做的是以下操作（使用伪代码）：


创建一个U大的矩阵size(Imlabel,1) x size(Imlabel,2) x num_class。该矩阵的每个切片将确定该切片中每个坐标(i,j)处的8个相邻像素的吉布斯能量成本。
对于每个类别标签x ...

一个。对于图像中的每个像素(i,j) ...


查找在(i,j)的8像素邻域中定义的等于x的那些位置，并将它们设置为beta
查找在(i,j)的8像素邻域中定义的不等于x的那些位置，并将它们设置为-beta
计算此8像素邻域内所有这些beta和-beta值的总和，并将其设置在位置U(i,j,x)



因此，我要做的是创建一个3D成本矩阵...让我们称之为C，其中每个切片y具有与size(Imlabel,1) x size(Imlabel,2)相同的尺寸，但是切片中的每个位置(i,j)将是 如果 beta，否则为 Imlabel(i,j) == y。完成此操作后，您可以对该矩阵执行N-D卷积，但是要使用2D内核，以便可以分别计算每个切片的8像素邻域之和。



当今的魔术菜单由 -beta， bsxfun和 convn的副食组成，以使锅变甜。

要获取成本矩阵 permute，您可以执行此操作。假设 C的边界没有填充零：

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;

Imlabel函数用于创建单个3D向量，其范围从1到您定义的多个类。之所以需要这样做，是因为 permute做所谓的广播。在给定矩阵 bsxfun的情况下会发生这种情况，使用3D向量将与 Imlabel或equals函数一起创建3D成本矩阵。此成本矩阵的每个部分都会为您提供与所讨论标签相等或不相等的位置。

您可以验证从 eq输出的locations矩阵的正确性：

>> 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

如您所见，对于每个切片，我们可以看到哪些位置共享该切片枚举的标签，以及哪些位置不共享同一标签。一旦有了这些位置，就需要将其强制转换为 bsxfun，以便可以将其修改为成本矩阵，其中 double是不等于特定切片标签的位置，而 -beta是。这是通过 beta调用后的两个赋值语句完成的。

有了这个成本矩阵后，就可以像以前一样填充边界，但是请像在代码中所做的那样，确保将边界设置为成本 bsxfun：

mask = zeros(3,3); mask(2,2) = 1;
Cpad = convn(C, mask);
Cpad(Cpad == 0) = -beta;

代码的第一行表示您在代码中定义的 -beta掩码。现在要独立填充每个切片以使其边界为零，我们可以使用 [0 0 0; 0 1 0; 0 0 0]来帮助我们做到这一点。



现在我们已经完成了设置，您要做的就是每个切片独立执行Gibbs能量的计算：

mask2 = ones(3,3); mask2(2,2) = 0;
out = convn(Cpad, mask2, 'valid');

代码的第一行表示一个掩码，其中每个值（中间值除外）均为1，即为0，并且是3 x 3掩码。您要执行此操作的原因是因为它替换了双 convn循环逻辑中的 if/else语句。实际上，您正在做的是一个卷积，其中将除中间中间之外的所有8个邻居相加。这可以通过使用除中心像素以外的全1的掩码并将其设置为0来实现。

接下来，我们使用 for独立地计算每个切片的吉布​​斯能量，但是使用 convn标志，以便我们可以删除开头包含的零填充边界。

现在，如果您看一下 'valid'，它将与您使用 out的情况进行比较，但是由于索引到 U的方式略有变化。不过，使用 U的上述输出，您可以验证我们所拥有的是否正确，并且可以很好地处理边界情况：

>> 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，因此自 out以来，我们的总和为3，但是然后还有其他5个我们未考虑的方向，在您的代码中，将这些方向设置为 beta = 1，因此将其设置为 -beta。

我们还要看一下第六行第四列，看一下第二个标签或切片。这意味着任何等于标签2的基本方向，我们都应将其与 3 - 5 = -2相加，而与标签2不相等的任何基本方向均应为 beta / 1。如您所见，正好有4个标签为1，4个标签为2，因此这些标签应抵消并给我们0。

您可以验证对矩阵中所有其他位置所做的操作是否得出正确的计算结果。

完整的代码很简单：

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');

在定时的情况下，我们可以检查一下上述方法与原始循环代码的比较速度。我使用 -beta / -1来帮助实现这一点。

这是我设置的脚本，该脚本设置了所有相关变量，并测量了您的代码和我的代码所花费的时间：

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

运行上面的脚本，我有时会得到这些：

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的因数，这是可怕的。但是，您显示的是这么小的数据集。我们应该尝试使用更大的数据集，看看是否仍然有效。让我们制作具有10个类别标签的2500 x 2500矩阵。我将用以下内容替换 timeit矩阵：

rng(123); %// Set seed for reproducibility
num_class = 10;
Imlabel_orig = randi(10,2500,2500);

这将在2500 x 2500矩阵中产生1到10的随机整数。这样做并再次运行代码将得到：

Loop code time: 17.553669 seconds
Vectorized code time: 3.321950 seconds
Speedup factor: 5.284146

是的...好多了。在10个标签的2500 x 2500矩阵上，计算大约需要17.5秒，而带有 Imlabel的矢量化代码的性能要比原始循环代码高出近5.2倍。



至于归纳到吉布斯能量计算中的归一化项，您目前具有以下条件：

P=zeros(size(U));
Z=sum(exp(-U),3);
for i=1:num_class
   P(:,:,i)=exp(-U(:,:,i))./Z;
end

如果利用 bsxfun / convn / permute的优势，则只需执行以下操作：

Z = sum(exp(-out),3)); %// out is U in my code
P = bsxfun(@rdivide, exp(-out), Z);

这可以实现代码所执行的相同操作。对于 bsxfun中的每个切片，我们找到 P并将该切片取反作为输入，然后将每个术语除以 exp。对于 Z，将复制 bsxfun矩阵的切片数与 Z中的切片数相同，并且按元素进行除法非常类似于循环代码。