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

来源:互联网 发布:梦幻西游网络出错 编辑:程序博客网 时间:2024/03/29 14:38

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

图像背景分离也称图像切割,原理是利用聚类算法将图像进行聚类,聚为不同的若干个类,这些类别即包含了图像的前景和背景。
此项目包含以下几个方面

聚类方法:实现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    endend


 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 deleted 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; %yend


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    endend


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时,从一副图像将前景的猫分离出来,拼接到另一幅背景为草地的图像上,效果如图:


 

0 0
原创粉丝点击