使用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入门程序介绍*/