CS131-PA2 通过聚类实现前/背景分离 Foreground-Background Segmentation via Clustering

通过聚类实现前/背景分离

图像背景分离也称图像切割，原理是利用聚类算法将图像进行聚类，聚为不同的若干个类，这些类别即包含了图像的前景和背景。

此项目包含以下几个方面

聚类方法：实现K-Means和Hierarchical Agglomerative算法

点特征向量：实现ComputePositionColorFeatures（颜色和位置）方法和特征归一化

不同参数实验并分析

 

1.实现K-Means和Hierarchical Agglomerative算法

1.1 K-Means 是一种无监督的聚类算法，用来对样本数据进行聚类运算，k为聚类的数量，means平均，不断取离种子点中心最近均值

     根据k，对样本集中随机抽取k个样本，作为种子点

     计算所有样本到这k个种子点的距离，离某个种子点得最近的样本即划归为此类

    根据产生的k个数据集，重新计算新的种子点，可以使用数据集中样本的各个维度的均值

    不断迭代，当此k个种子点不再变化（稳定）或超过迭代次数阀值的时候，停止算法

一般用于数据挖掘，文档归类，未知类别数量的时候使用

详见百科：

http://baike.baidu.com/link?url=KAJX8qy8GtKzLiriueWdrr76LWbgVPPxDqcwPzaMBupH7-QJGZ34lykeISFY7EBBOisDIeH-BqMihqhvcWz6aK

具体实现：

function idx = KMeansClustering(X, k, visualize2D, centers)
% Run the k-means clustering algorithm.
%
% INPUTS
% X - An array of size m x n containing the points to cluster. Each row is
%     an n-dimensional point, so X(i, :) gives the coordinates of the ith
%     point.
% k - The number of clusters to compute.
% visualize2D - OPTIONAL parameter telling whether or not to visualize the
%               progress of the algorithm for 2D data. If not set then 2D
%               data will not be visualized.
% centers - OPTIONAL parameter giving initial centers for the clusters.
%           If provided, centers should be a k x n matrix where
%           centers(c, :) is the center of the cth cluster. If not provided
%           then cluster centers will be initialized by selecting random
%           rows of X. You don't need to use this parameter; it is mainly
%           here to make your code more easily testable.
%
% OUTPUTS
% idx     - The assignments of points to clusters. idx(i) = c means that the
%           point X(i, :) has been assigned to cluster c.

if ~isa(X, 'double')
X = double(X);
end
m = size(X, 1);
n = size(X, 2);

% If we are going to display the clusters graphically then create a
% figure in which to draw the visualization.
figHandle = [];
if n == 2 && visualize2D
figHandle = figure;
end

% If initial cluster centers were not provided then initialize cluster
% centers to random rows of X. Each row of the centers variable should
% contain the center of a cluster, so that centers(c, :) is the center
% of the cth cluster.
if ~exist('centers', 'var')
centers = zeros(k, n);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
random_row = randperm(m);%随机1到m的整数
random_row = random_row(1:k);%取前k个数
centers=X(random_row,:);%随机赋给中心点
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
end

% The assignments of points to clusters. If idx(i) == c then the point
% X(i, :) belongs to the cth cluster.
idx = zeros(m, 1);

% The number of iterations that we have performed.
iter = 0;

% If the assignments of points to clusters have not converged after
% performing MAX_ITER iterations then we will break and just return the
% current cluster assignments.
MAX_ITER = 100;

while true
% Store old cluster assignments
old_idx = idx;

% Compute distances from each point to the centers and assign each
% point to the closest cluster.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for i = 1 : m
dists = [];
for j = 1:k
dist = sum((X(i,:) - centers(j,:)).^2); %循环所有样本点，每个样本点计算和中心点的距离
dists = [dists,dist];%保存和中心点的距离
end
[min_dist,idx(i)] = min(dists); %取最小距离的中心作为同一类
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Break if cluster assignments didn't change
if idx == old_idx
break;
end

% Update the cluster centers
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%更新新的中心点，取各类x,y的均值作为新的中心点
for i = 1:k
cluster_i = X(find(idx==i),:);%i聚类的所有样本点
centers(i,:) = mean(cluster_i); %取各类x,y的均值作为新的中心点
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Display the points in the 2D case.
if n == 2 && visualize2D
VisualizeClusters2D(X, idx, centers, figHandle);
pause;
end

% Stop early if we have performed more than MAX_ITER iterations
iter = iter + 1;
if iter > MAX_ITER
break;
end
end
end

 1.2 Hierarchical Agglomerative算法 层次聚类，

原理：设定一个最终要聚成的类别数k，初始化时每个样本点为一个类，计算每个样本对应最小距离的样本（类），合并这两个类，直到类别数小于或等于k

实现如下：

function idx = HAClustering(X, k, visualize2D)
% Run the hierarchical agglomerative clustering algorithm.
%
% The algorithm is conceptually simple:
%
% Assign each point to its own cluster
% While the number of clusters is greater than k:
%   Compute the distance between all pairs of clusters
%   Merge the pair of clusters that are closest to each other
%
% There are many valid ways of determining the distance between two
% clusters. For this assignment we will define the distance between two
% clusters to be the Euclidean distance between their centroids.
%
% Recomputing the centroids of all clusters and the distances between all
% pairs of centroids at each step of the loop would be very slow.
% Thankfully most of the distances and centroids remain the same in
% successive iterations of the outer loop; therefore we can speed up the
% computation by only recomputing the centroid and distances for the new
% merged cluster.
%
% Even with this trick, this algorithm will consume a lot of memory and run
% very slowly when clustering large sets of points. In practice, you
% probably do not want to use this algorithm to cluster more than 10,000
% points.
%
% NOTE: To get full credit for this part you should only update the cluster
% centroids and distances for the merged and delet
e9c7
ed clusters at each
% iteration of the main loop.
%
% INPUTS
% X - An array of size m x n containing the points to cluster. Each row is
%     an n-dimensional point, so X(i, :) gives the coordinates of the ith
%     point.
% k - The number of clusters to compute.
% visualize2D - OPTIONAL parameter telling whether or not to visualize the
%               progress of the algorithm for 2D data. If not set then 2D
%               data will not be visualized.
%
% OUTPUTS
% idx     - The assignments of points to clusters. idx(i) = c means that the
%           point X(i, :) has been assigned to cluster c.

if nargin < 3
visualize2D = false;
end

if ~isa(X, 'double')
X = double(X);
end

% The number of points to cluster.
m = size(X, 1);

% The dimension of each point.
n = size(X, 2);

% The number of clusters that we currently have.
num_clusters = m;

% The assignment of points to clusters. If idx(i) = c then X(i, :) is
% assigned to cluster c. Since each point is initally assigned to its
% own cluster, we initialize idx to the column vector [1; 2; ...; m].
idx = (1:m)';

% The centroids of all clusters. The row centroids(c, :) is the
% centroid of the cth cluster. Since each point starts in its own
% cluster, the centroids are initially the same as the data matrix.
centroids = X;

% The number of points in each cluster. If cluster_sizes(c) = n then
% cluster c contains n points. Since each point is initially assigned
% to its own cluster, each cluster size is initialized to 1.
cluster_sizes = ones(m, 1);

% The distances between the centroids of the clusters. For ci != cj,
% dists(ci, cj) = d means that the Euclidean distance between
% centroids(ci, :) and centroids(cj, :) is d. We will choose the pair
% of clusters to merge at each step by finding the smallest element of
% the dists matrix. We never want to merge a cluster with itself;
% therefore we set the diagonal elements of dists to be +Inf.
dists = squareform(pdist(centroids)); %pdist是求矩阵内行与行之间的欧式距离,返回1行向量，使用squareform转换为原距离（centroids维度一样）的矩阵
dists = dists + diag(Inf(m, 1));%将dists矩阵对角置1，因为不需要跟自身合并

% If we are going to display the clusters graphically then create a
% figure in which to draw the visualization.
figHandle = [];
if n == 2 && visualize2D
figHandle = figure;
end

% In the 2D case we can easily visualize the starting points.
if n == 2 && visualize2D
VisualizeClusters2D(X, idx, centroids, figHandle);
pause;
end

while num_clusters > k

% Find the pair of clusters that are closest together.
% Set i and j to be the indices of the nearest pair of clusters.
i = 0;
j = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%每次都找最近的两个类？
[min_dists,columns]=min(dists); %最小值所在的列索引和值
[min_dist,row] = min(min_dists); %从最小值所在的列索引和值中获取最小行索引
i = row;
j = columns(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Make sure that i < j
if i > j
t = i;
i = j;
j = t;
end

% Next we need to merge cluster j into cluster i.
%
% We also need to 'delete' cluster j. We will do this by setting
% its cluster size to 0 and moving its centroid to +Inf. This
% ensures that the distance from cluster j to any other cluster is
% +Inf.

% Assign all points currently in cluster j to cluster i.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%分配（指定）j类的点到i类，为了合并做准备
idx(find(idx==j)) = idx(i);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Compute the new centroid for clusters i and set the centroid of
% cluster j to +Inf.
% HINT: You should be able to compute both updated cluster
% centroids in O(1) time.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
centroid_i = centroids(i,:); %当前i类质心
centroid_j = centroids(j,:); %当前j类质心
centroid_ij = [centroid_i;centroid_j];
centroid_new = mean(centroid_ij); %平均值求新的质心
centroids(i,:) = centroid_new; %平均值求新的质心
centroids(j,:) = Inf;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Update the size of clusters i and j.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
cluster_sizes(i)=cluster_sizes(i)+cluster_sizes(j);
cluster_sizes(j) = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% Update the dists array. In particular, we need to update the
% distances from clusters i and j to all other clusters.
% Hint: You might find the pdist2 function useful.
% Hint: Remember that the diagonal of dists must be +Inf.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            YOUR CODE HERE                           %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%i和j类合并为i类之后，i类到其他类的距离
%Group Average：这种方法看起来相对有道理一些，也就是把两个集合中的点两两的距离全部放在一起求一个平均值，相对也能得到合适一点的结果。
%参考http://blog.pluskid.org/?p=407

dists(i,:) = (dists(i,:)+dists(j,:))/2;
%dists(i,:) = pdist2(centroids(i,:),centroids(:,:)); %新的质心到其他类质心的距离
dists(i,i) = Inf;
dists(j,:) = Inf; %j类到所有类的距离清空
dists(:,j) = Inf; %所有类到j类的距离清空
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                     %
%                            END YOUR CODE                            %
%                                                                     %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% If everything worked correctly then we have one less cluster.
num_clusters = num_clusters - 1;

% In the 2D case we can easily visualize the clusters.
if n == 2 && visualize2D
VisualizeClusters2D(X, idx, centroids, figHandle);
pause;
end

end

% Now we need to reindex the clusters from 1 to k
idx = ReindexClusters(idx);
end

2.点特征向量

众所周知，每幅图像都是由像素点组成的。想象一下，怎么区分像素和像素是属于同一类的呢？按照常规思维，我们会把颜色相同，或者位置相近的像素点视为同一类。

2.1 Color Features 颜色特征，这个很简单，就是直接每个像素取RGB3个维度的值作为特征向量

2.2Color and Position Features 颜色结合位置特征，在Color Features的基础上增加位置信息

代码：

function features = ComputePositionColorFeatures(img)
% Compute a feature vector of colors and positions for all pixels in the
% image. For each pixel in the image we compute a feature vector
% (r, g, b, x, y) where (r, g, b) is the color of the pixel and (x, y) is
% its position within the image.
%
% INPUT
% img - Array of image data of size h x w x 3.
%
% OUTPUT
% features - Array of computed features of size h x w x 5 where
%            features(i, j, :) is the feature vector for the pixel
%            img(i, j, :).

height = size(img, 1);
width = size(img, 2);
features = zeros(height, width, 5);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                              YOUR CODE HERE                             %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
features(:,:,1) = img(:,:,1); %r
features(:,:,2) = img(:,:,2); %g
features(:,:,3) = img(:,:,3); %b

%img_x = zeros(height,width);
index_x = [1:height]; %x坐标
index_y = [1:width]; %y坐标
img_x = (ones(width,1)*[1:height])';
img_y = ones(height,1)*[1:width];

features(:,:,4) = img_y; %x
features(:,:,5) = img_x; %y
end

2.3 特征归一化

零均值单位方差zero mean and unit variance

均值，顾名思义一向量的平均值，matlab中，若是m*n矩阵，则使用mean函数，可以直接计算

方差，定义：u^2 = （(x1-x)^2+(x2-x)^2+....+(xn-x)^2）/n-1 其中x为均值，意义是表示数据的离散程度

距离，如两组数据 甲100万，乙20万，平均60万，丙70万，丁50万，平均也是60万

但是甲乙的方差比丙丁的方差大，说明波动较大，也就是越不“平均”

zero mean and unit variance 公式

feature = (fi-x)/u

由于特征各个维度数据范围不同，来源不同，所以为了更好的进行数据分析，所以对特征数量进行标准化

实现代码：

function featuresNorm = NormalizeFeatures(features)
% Normalize image features to have zero mean and unit variance. This
% normalization can cause k-means clustering to perform better.
%
% INPUTS
% features - An array of features for an image. features(i, j, :) is the
%            feature vector for the pixel img(i, j, :) of the original
%            image.
%
% OUTPUTS
% featuresNorm - An array of the same shape as features where each feature
%                has been normalized to have zero mean and unit variance.

features = double(features);
featuresNorm = features;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                              %
%                                YOUR CODE HERE:                               %
%                                                                              %
%                HINT: The functions mean and std may be useful                %
%                                                                              %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[height,width,dim]=size(featuresNorm);
sum_u = zeros(1,dim);
u = zeros(1,dim);
%    %传统循环计算均值
%     for i = 1:height
%         for j = 1:width
%             sum_u = sum_u+reshape(features(i,j,:),1,dim);
%         end
%     end
%     u = sum_u/(height*width); %均值

%使用mean函数
feat = reshape(featuresNorm,height*width,dim); %将100*200*10的矩阵转换为20000*10的矩阵，方便使用mean
u = mean(feat); %返回一个row_vector 每个值为列值的平均

%     %传统循环计算方差
%     sum_v = zeros(1,dim);
%     v = zeros(1,dim);
%     for i = 1:height
%         for j = 1:width
%             feat = features(i,j,:);
%             feat = reshape(feat,1,dim);
%             sum_v = sum_v+(feat-u).^2;
%         end
%     end
%     v = (sum_v/(height*width-1)); %方差
%     v = v.^(1/2);
v = std(feat); %直接调用std函数计算方差

for i = 1:height
for j = 1:width
featuresNorm(i,j,:) = (reshape(featuresNorm(i,j,:),1,dim)-u)./v;
end
end
end

2.4 其他特征

可以尝试计算多种不同的特征，并结合起来行程新的特征，用于提高聚类的效果。这里额外计算了梯度和赋值结合起来的特征：

function features = ComputeFeatures(img)
% Compute a feature vector for all pixels of an image. You can use this
% function as a starting point to implement your own custom feature
% vectors.
%
% INPUT
% img - Array of image data of size h x w x 3.
%
% OUTPUT
% features - Array of computed features for all pixels of size h x w x f
%            such that features(i, j, :) is the feature vector (of
%            dimension f) for the pixel img(i, j, :).

img = double(img);
height = size(img, 1);
width = size(img, 2);
features = zeros([height, width, 1]);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                                                                         %
%                              YOUR CODE HERE                             %
%                                                                         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
features = zeros([height, width, 5]); %重新定义为5个维度
gray_img = rgb2gray(img); %转化为灰度图像
%颜色
features(:,:,1) = img(:,:,1); %r
features(:,:,2) = img(:,:,2); %g
features(:,:,3) = img(:,:,3); %b

%计算梯度
dx = [-1 0 1];%x方向滤窗
dy = [-1;0;1];%y方向滤窗
img_dx = filter2(dx,gray_img);%x方向做点积,默认参数是same
img_dy = filter2(dy,gray_img);%y方向做点积

%计算梯度的幅值和幅角
grad_mag = (img_dx.^2+img_dy.^2).^(1/2);
grad_theta = atan2(img_dy,img_dx);

features(:,:,4) = grad_mag; %幅值
features(:,:,5) = grad_theta; %幅角
end

 

3.实验结果

初始代码已经实现了分离，我们只需提供聚类方法、类数、特征类型等参数即可看到实验结果

当k=10，使用k-means聚类方法，归一化为true 使用color_feature和color_position_feature对应的结果为：

color_feature:



 

color_position_feature:



说明在颜色特征基础上加了位置信息后，聚类的准确率得到了提升

 

最后，使用k=4,特征为color_feature,聚类方法为k-means时，从一副图像将前景的猫分离出来，拼接到另一幅背景为草地的图像上，效果如图：