使用k-means对3D网格模型进行分割

由于一些原因，最近在做网格分割的相关工作。网格分割的方法有很多，如Easy mesh cutting、K-means、谱分割、基于SDF的分割等。根据对分割要求的不同，选取合适的分割方法。本文中使用了较为简单的k-means对网格进行分割。

K-means原理

K-means是一种简单的聚类方法，聚类属于无监督学习，聚类的样本中却没有给定y，只有特征x，比如假设宇宙中的星星可以表示成三维空间中的点集(x,y,z)。聚类的目的是找到每个样本x潜在的类别y，并将同类别y的样本x放在一起。对于上述的星星，聚类后结果是一个个星团，星团里面的点相互距离比较近，不同星团间的星星距离就比较远了。

算法描述

（１）从数据集中随机抽取k个质心作为初始聚类的中心；
（２）计算数据集中所有的点到这k个点的距离，将点归到离其最近的聚类里；
（３）调整聚类中心，即将聚类的中心移动到聚类的几何中心（即平均值）处；
（４）重复第2步和第3步，直到聚类的中心不再移动，此时算法收敛。

matlab代码如下

function cluster_labels = k_means(data, centers, num_clusters)

%K_MEANS Euclidean k-means clustering algorithm.

%

%   Input    : data           : N-by-D data matrix, where N is the number of data,

%                               D is the number of dimensions

%              centers        : K-by-D matrix, where K is num_clusters, or

%                               'random', random initialization, or

%                               [], empty matrix, orthogonal initialization

%              num_clusters   : Number of clusters

%

%   Output   : cluster_labels : N-by-1 vector of cluster assignment

% Parameter setting

%

iter = 0;

qold = inf;

threshold = 0.001;

%

% Check if with initial centers

%

if strcmp(centers, 'random')

  disp('Random initialization...');

  centers = random_init(data, num_clusters);

elseif isempty(centers)

  disp('Orthogonal initialization...');

  centers = orth_init(data, num_clusters);

end

%

% Double type is required for sparse matrix multiply

%

data = double(data);

centers = double(centers);

%

% Calculate the distance (square) between data and centers

%

n = size(data, 1);

x = sum(data.*data, 2)';

X = x(ones(num_clusters, 1), :);

y = sum(centers.*centers, 2);

Y = y(:, ones(n, 1));

P = X + Y - 2*centers*data';

%

% Main program

%

while 1

  iter = iter + 1;

  % Find the closest cluster for each data point

  [val, ind] = min(P, [], 1);

  % Sum up data points within each cluster

  P = sparse(ind, 1:n, 1, num_clusters, n);

  centers = P*data;

  % Size of each cluster, for cluster whose size is 0 we keep it empty

  cluster_size = P*ones(n, 1);

  % For empty clusters, initialize again

  zero_cluster = find(cluster_size==0);

  if length(zero_cluster) > 0

    disp('Zero centroid. Initialize again...');

    centers(zero_cluster, :)= random_init(data, length(zero_cluster));

    cluster_size(zero_cluster) = 1;

  end

  % Update centers

  centers = spdiags(1./cluster_size, 0, num_clusters, num_clusters)*centers;

  % Update distance (square) to new centers

  y = sum(centers.*centers, 2);

  Y = y(:, ones(n, 1));

  P = X + Y - 2*centers*data';

  % Calculate objective function value

  qnew = sum(sum(sparse(ind, 1:n, 1, size(P, 1), size(P, 2)).*P));

  mesg = sprintf('Iteration %d:\n\tQold=%g\t\tQnew=%g', iter, full(qold), full(qnew));

  disp(mesg);

  % Check if objective function value is less than/equal to threshold

  if threshold >= abs((qnew-qold)/qold)

    mesg = sprintf('\nkmeans converged!');

    disp(mesg);

    break;

  end

  qold = qnew;

end

cluster_labels = ind';

%-----------------------------------------------------------------------------

function init_centers = random_init(data, num_clusters)

%RANDOM_INIT Initialize centroids choosing num_clusters rows of data at random

%

%   Input : data         : N-by-D data matrix, where N is the number of data,

%                          D is the number of dimensions

%           num_clusters : Number of clusters

%

%   Output: init_centers : K-by-D matrix, where K is num_clusters

rand('twister', sum(100*clock));

init_centers = data(ceil(size(data, 1)*rand(1, num_clusters)), :);

function init_centers = orth_init(data, num_clusters)

%ORTH_INIT Initialize orthogonal centers for k-means clustering algorithm.

%

%   Input : data         : N-by-D data matrix, where N is the number of data,

%                          D is the number of dimensions

%           num_clusters : Number of clusters

%

%   Output: init_centers : K-by-D matrix, where K is num_clusters

%

% Find the num_clusters centers which are orthogonal to each other

%

Uniq = unique(data, 'rows'); % Avoid duplicate centers

num = size(Uniq, 1);

first = ceil(rand(1)*num); % Randomly select the first center

init_centers = zeros(num_clusters, size(data, 2)); % Storage for centers

init_centers(1, :) = Uniq(first, :);

Uniq(first, :) = [];

c = zeros(num-1, 1); % Accumalated orthogonal values to existing centers for non-centers

% Find the rest num_clusters-1 centers

for j = 2:num_clusters

  c = c + abs(Uniq*init_centers(j-1, :)');

  [minimum, i] = min(c); % Select the most orthogonal one as next center

  init_centers(j, :) = Uniq(i, :);

  Uniq(i, :) = [];

  c(i) = [];

end

clear c Uniq;

网格分割

对网格使用K-means进行聚类，首先要构造用于聚类的特征，也就是要构造data、center和num_clusters
如果只使用网格顶点坐标作为聚类特征，并使用'random'初始化聚类中心，聚类结果如下：

模型顶点：53054，面：106104，聚类特征：顶点坐标，聚类数目：4
       

如果将顶点法向也加入聚类的特征中，得到的结果如下：
      

附获取顶点法向的代码如下：

/*获取所有顶点的法向*/

    _mesh->request_face_normals();

    _mesh->request_vertex_normals();

    _mesh->update_normals();

    int vnidx = 0;

    for (auto vit = _mesh->vertices_begin(); vit != _mesh->vertices_end(); ++vit,vnidx++)

    {

        auto vertex = vit.handle();

        OpenMesh::Vec3d v_normal;

         v_normal = _mesh->normal(vertex);

        Vnx[vnidx] = v_normal.data()[0];

        Vny[vnidx] = v_normal.data()[1];

        Vnz[vnidx] = v_normal.data()[2];

    }

    /*一部分参考YQ，一部分参考OpenMesh入门程序介绍*/