使用高斯混合模型对图像进行聚类

Posted 2023-03-12

【中文标题】使用高斯混合模型对图像进行聚类

【英文标题】：Clustering an image using Gaussian mixture models

【发布时间】：2015-05-09 12:55:19

【问题描述】：

我想使用 GMM（高斯混合模型对二值图像进行聚类，并且还想在二值图像本身上绘制聚类质心。

我以此为参考： http://in.mathworks.com/help/stats/gaussian-mixture-models.html

这是我的初始代码

 I=im2double(imread('sil10001.pbm'));
  K = I(:);
  mu=mean(K);
  sigma=std(K);
  P=normpdf(K, mu, sigma);
   Z = norminv(P,mu,sigma);
  X = mvnrnd(mu,sigma,1110);
  X=reshape(X,111,10);


 scatter(X(:,1),X(:,2),10,'ko');

options = statset('Display','final');
gm = fitgmdist(X,2,'Options',options);



idx = cluster(gm,X);
cluster1 = (idx == 1);
cluster2 = (idx == 2);


 scatter(X(cluster1,1),X(cluster1,2),10,'r+');
 hold on

  scatter(X(cluster2,1),X(cluster2,2),10,'bo');
  hold off
  legend('Cluster 1','Cluster 2','Location','NW')


  P = posterior(gm,X);

 scatter(X(cluster1,1),X(cluster1,2),10,P(cluster1,1),'+')
 hold on
 scatter(X(cluster2,1),X(cluster2,2),10,P(cluster2,1),'o')
 hold off
 legend('Cluster 1','Cluster 2','Location','NW')
 clrmap = jet(80); colormap(clrmap(9:72,:))
 ylabel(colorbar,'Component 1 Posterior Probability')

但问题是我无法在主二进制图像中绘制从 GMM 接收到的簇质心。我该怎么做？

**现在假设我在一个序列中有 10 个这样的图像并且我想将它们的平均位置信息存储在两个单元格数组中，那么我该怎么做。这是我的新问题的代码**

    images=load('gait2go.mat');%load the matrix file
    for i=1:10

   Ii=images.resulti;
  Ii=im2double(Ii);

   %determine 'white' pixels, size of image can be [M N], [M N 3] or [M N 4]
  Idims=size(Ii);
  whites=true(Idims(1),Idims(2));

    df=Ii;
      %we add up the various color channels
 for colori=1:size(df,3)
  whites=whites & df(:,:,colori)>0.5;
 end

%choose indices of 'white' pixels as coordinates of data
[datax datay]=find(whites);

%cluster data into 10 clumps
  K = 10;               % number of mixtures/clusters
  cInd = kmeans([datax datay], K, 'EmptyAction','singleton',...
   'maxiter',1000,'start','cluster');

%get clusterwise means
 meanx=zeros(K,1);
 meany=zeros(K,1);  
  for i=1:K
   meanx(i)=mean(datax(cInd==i));
   meany(i)=mean(datay(cInd==i));

 end

 xci=meanx(i);%cell array contaning the position of the mean for the 10    
 images
  xbi=meany(i);

figure;
gscatter(datay,-datax,cInd); %funky coordinates for plotting according to      
 image
 axis equal;
  hold on;
  scatter(meany,-meanx,20,'+'); %same funky coordinates


 end

我能够分割出 10 张图像，但没有存储在单元格数组 xc 和 xb 中的均值值。它们仅存储 [] 代替均值值

【问题讨论】：

看来您设置了一个很好的算法。您在将此算法转换为代码时遇到问题吗？到目前为止，您尝试过什么吗？你应该知道 *** 不是一个代码编写服务。您需要向我们展示您迄今为止所做的尝试，我们会帮助您解决问题。毫不费力的问题将不可避免地被关闭。

@user3371423 请格式化您的代码。

那些圆圈不是散点图的结果吗？我认为您只是将所有颜色更改为黑白，这就是为什么您什么都看不到。删除 clrmap 行

@AnderBiguri 这是一个不同代码的输出，它是 k 意味着不是这个我想要使用 GMM 的输出

【参考方案1】：

我决定发布您的问题的答案（您的问题由最大似然猜测确定：P），但我写了一个详尽的介绍。请仔细阅读，因为我认为你很难理解你想使用的方法，你也很难理解为什么别人不能用你通常的提问方式来帮助你。您的问题有几个问题，包括与代码相关的和概念性的。让我们从后者开始。

有问题的问题

您说您想使用高斯混合建模对图像进行聚类。虽然我通常不熟悉集群，但在查看了您的 reference 和 the wonderful SO answer you cited elsewhere（以及来自 @rayryeng 的快速 101）之后，我认为您完全走错了路。

顾名思义，高斯混合建模使用混合高斯（即正态）分布对您的数据集进行建模。这种方法流行的原因是，当你测量各种量的时候，很多时候你会发现你的数据大部分是像正态分布一样分布的（这实际上是它被称为normal)。这背后的原因是central limit theorem，这意味着在许多情况下，合理独立的随机变量的总和往往是正常的。

现在，聚类，另一方面，只是意味着根据某些标准将您的数据集分成不相交的更小的束。主要标准通常是（某种）距离，因此您希望在较大的数据集中找到“紧密的数据块”。您通常需要在执行 GMM 之前对数据进行聚类，因为在不必猜测聚类的情况下找到数据背后的高斯已经足够困难了。我对所涉及的过程还不够熟悉，无法告诉 GMM 算法在你的原始数据上工作的效果如何（但我希望无论如何，许多实现都是从聚类步骤开始的）。

为了更接近您的问题：我猜您想做某种图像识别。查看图片，您想要获得更强烈相关的肿块。这是集群。如果你看一张动物园的照片，你会看到，比如说，一头大象和一条蛇。两者都有其独特的形状，并且彼此很好地分开。如果您对图像进行聚类（并且蛇没有骑大象，neither did it eat it），您会发现两个块：一个是大象形的，一个是蛇形的。现在，在这些数据集上使用 GMM 是没有意义的：大象，尤其是蛇，不像多元高斯分布。但是，如果您只想知道不同动物在您的照片中的位置，您一开始就不需要这个。

仍然使用示例，您应该确保将数据聚集到适当数量的子集中。如果您尝试将您的动物园图片分成 3 个集群，您可能会得到第二条假蛇：大象的鼻子。随着集群数量的增加，您的分区可能越来越没有意义。

你的方法

你的代码没有给你任何合理的东西，这有一个很好的理由：它从一开始就没有意义。看开头：

I=im2double(imread('sil10001.pbm'));
K = I(:);
mu=mean(K);
sigma=std(K);
X = mvnrnd(mu,sigma,1110);
X=reshape(X,111,10);

您读取二进制图像，将其转换为双精度图像，然后将其拉伸成向量并计算该向量的均值和偏差。您基本上将您的内部图像涂抹成 2 个值：平均强度和偏差。然后你用这些参数生成111*10标准法线点，并尝试在前两组111上做GMM。它们都是具有相同参数的独立法线。因此，您可能会得到两个重叠的高斯分布，均值相同且偏差相同。

我认为您在网上找到的示例让您感到困惑。当你做 GMM 时，你已经有了你的数据，所以不应该涉及伪正态数。但是当人们发布示例时，他们也会尝试提供可重复的输入（嗯，他们中的一些人会这样做，nudge nudge wink wink）。一个简单的方法是生成简单高斯的并集，然后可以将其输入 GMM。

所以，我的意思是，您不必生成随机数，而必须使用图像数据本身作为程序的输入。 而且您可能只想对图像进行聚类，而不是实际使用 GMM 在聚类上绘制土豆，因为您希望将身体部位聚类到关于人类的图像中。大多数身体部位的形状不像多元高斯分布（男性和女性有一些明显的例外）。

我认为你应该怎么做

如果您真的想对图像进行聚类，就像您添加到问题中的图中那样，那么您应该使用像 k-means 这样的方法。但是话又说回来，你已经有一个程序可以做到这一点，不是吗？所以我真的不认为我可以回答“如何将我的图像与 GMM 聚类？”的问题。相反，这是对“如何对图像进行聚类？”的答案。用k-means，但至少这里会有一段代码。

%set infile to what your image file will be
infile='sil10001.pbm';

%read file
I=im2double(imread(infile));

%determine 'white' pixels, size of image can be [M N], [M N 3] or [M N 4]
Idims=size(I);
whites=true(Idims(1),Idims(2));

%we add up the various color channels
for colori=1:Idims(3)
    whites=whites & I(:,:,colori)>0.5;
end

%choose indices of 'white' pixels as coordinates of data
[datax datay]=find(whites);

%cluster data into 10 clumps
K = 10;               % number of mixtures/clusters
cInd = kmeans([datax datay], K, 'EmptyAction','singleton',...
    'maxiter',1000,'start','cluster');

%get clusterwise means
meanx=zeros(K,1);
meany=zeros(K,1);
for i=1:K
    meanx(i)=mean(datax(cInd==i));
    meany(i)=mean(datay(cInd==i));
end

figure;
gscatter(datay,-datax,cInd); %funky coordinates for plotting according to image
axis equal;
hold on;
scatter(meany,-meanx,20,'ko'); %same funky coordinates

这就是它的作用。它首先像你一样读取你的图像。然后它尝试通过检查每个颜色通道（可以是 1、3 或 4）是否比 0.5 亮来确定“白色”像素。然后，您指向聚类的输入数据将是您的白色像素的x 和y“坐标”（即索引）。

接下来它通过kmeans 进行聚类。这部分代码松散地基于the already cited answer of Amro。我必须设置一个很大的最大迭代次数，因为从图片中没有 10 个清晰的簇的意义上说，这个问题是不适定的。然后我们计算每个集群的mean，并用gscatter绘制集群，用scatter绘制平均值。请注意，为了使图片在scatter 图中朝向正确的方向，您必须围绕输入坐标移动。或者，您可以在开头相应地定义datax 和datay。

这是我的输出，使用您在问题中提供的已处理数字运行：

【参考方案2】：

我相信你一定在情节中犯了一个幼稚的错误，这就是为什么你只看到一条直线：你只绘制了 x 值。

在我看来，scatter 命令中的第二个参数应该是X(cluster1,2) 或X(cluster2,2)，具体取决于代码中使用的是哪个scatter 命令。

【讨论】：

会起作用吗？因为 X 只有一列矩阵。这就是为什么在分散命令中我写 x(:1),x(:1) 而不是 x(:,1),x(:,2)

【参考方案3】：

代码可以更简单：

%read file

I=im2double(imread('sil10340.pbm'));
%choose indices of 'white' pixels as coordinates of data
[datax datay]=find(I);
%cluster data into 10 clumps
 K = 10;               % number of mixtures/clusters
[cInd, c] = kmeans([datax datay], K, 'EmptyAction','singleton',...
'maxiter',1000,'start','cluster');
 figure;
gscatter(datay,-datax,cInd); %funky coordinates for plotting according to    
image
axis equal;
hold on;
 scatter(c(:,2),-c(:,1),20,'ko'); %same funky coordinates

我认为不需要循环，因为 c 本身返回一个 10x2 双数组，其中包含方法的位置