点云-均值漂移-均适用

%实现均值漂移算法 function [idx, cluCenter, clucell] = meanshift(point, r, kernel) %input： %point：是输入的点 %r：是中心点扫描的半径，也可以叫做带宽 %kernel：是偏移值的计算方法，分为平均偏移值和核函数，核函数就是高斯分布 %output： % idx为每个点对应的类别 % clucount：每个类别包含的点的个数 % center：聚类中心点 %----------------------检查输入项-------------------- if nargin < 2 print('输入参数少于2个，至少包含点数据和所搜半径') end %----------------------初始化数据--------------------- numcluster = 0; %此行记录类别的数量 [n,m] = size(point); %记录点的行列 numn = n; %存储n的个数，其实也是指针 initpind = 1:n; %点的序列 bandsq = r^2; %计算带宽的平方,后面计算的距离值也是平方，所以正好对应 stopthr = 1e-3*r; %建立漂移的阈值，当大于阈值时停止 cluCenter = []; %建立存储中心点的列表 beenvisited = false(n,1); %建立一列表示是否被访问过，0表示没有访问过，1表示访问过 cluvotes = zeros(n,1,'uint16'); %建立类别投票项 clusMemcell = []; %建立存储元组存放每个类点的索引 %---------------------选择偏移函数-------------------- switch kernel case 'flat' %使用平均函数 kerfun = @(x,d) mean(x,1); case 'gaussian' %使用高斯核函数 kerfun = @(x,d) gaussm(x,d,r); otherwise %可以在此处扩展新的核函数 error("选择flat或是gaussian"); end %--------------------MS的核心代码--------------------------- %---------------首先先挑选种子点---------------------------- while numn %区域生长分割这里使用for循环，这里使用了while循环。 %但是其实两者都一样，前者有固定点数循环 %后者循环到所有点的被访问一遍 tempind = ceil((numn-1e-6)*rand); %这里n表示点的个数，即点的行数，rand是随机0—1之间的小数 %通过rand*点的个数，随机获取种子点的索引 %在原来的代码中使用了ceil，表示向上取整 %本代码中使用round四舍五入 seedind = initpind(tempind); %根据一开始做点数量列表，得到种子点的索引 mymean = point(seedind,:); %根据索引值得到种子点 mymembers = []; %将归属类别的点进行存储，mymember就是这样的空列表 thisvotes = zeros(n,1,'uint16'); %建立临时投票列 %-------------开始聚类---------------- while true %从这里就是开始聚类了 sqdisall = sum(bsxfun(@minus,mymean,point).^2,2); %计算种子点到其他点的距离的平方 %sum(,2)中2表示按行计算和 %bsxfun是两个功能的函数，其中minus表示mymean-point %由于高斯核函数需要用到距离，所以不能用临近搜索 inidx = find(sqdisall < bandsq); %这里寻找在带宽r内的点，返回索引 thisvotes(inidx) = thisvotes(inidx)+1; %这里建立投票，为带宽r内的点标记投票，标记类别 myoldmean = mymean; %存储种子点，保存为老的种子点 mymean = kerfun(point(inidx,:),sqrt(sqdisall(inidx))); %sqrt求的就是距离的开方，这一步是求的新种子点 mymembers = [mymembers;inidx]; %存储了类的点 %这一步的目的是储存初始类别的索引，在后面的步骤进行迭代 %收敛。不直接用索引的原因是每次索引到的点不一致，在第一次索引后， %后面搜索到的索引数量不会改变，会影响聚类。在每次聚类后， %都会回到mymembers=[]，可以往复储存 beenvisited(mymembers) = true; %将点标记为已访问 %-------------判断聚类是否收敛----------------- if norm(mymean - myoldmean) < stopthr %这里判断新种子点移动的距离是否在阈值之内，不满足阈值 mergew = 0; %创建类别、种子的指针 for cn = 1:numcluster %第一次聚类中numcluster的值为0，所以是1：0，不循环 disother = norm(mymean-cluCenter(cn,:)); %循环计算种子点到其他已存在种子点的长度，距离 if disother < r/2 mergew = cn; break; %这一步主要是通过种子点和已存在的种子点做对比 %判断两者是否相近，如果相近，则可以认为是一个类别 %计算两个相近种子点的权值，更新种子点 %mergew表示了种子点的索引 end end %----------更新种子点------------- %-----------计算权重-------------- if mergew > 0 %首次聚类时只有一个种子点，所有该判断不生效 nc = numel(mymembers); %计算当前种子点包含类别中点的数量 no = numel(clusMemcell{mergew}); %这里是统计相似种子点的类别中的点的数量 nw = [nc;no]/(nc+no); %得到对应点的权重信息 clusMemcell{mergew} = unique([clusMemcell{mergew};mymembers]); %更新相似种子点的成员，这里面的值是索引号 cluCenter(mergew,:) = nw(1)*mymean + nw(2)*myoldmean; %更新相似种子点的坐标，根据权重更新 cluvotes(:,mergew) = cluvotes(:,mergew) + thisvotes; %更新投票项 else numcluster = numcluster+1; %类别+1 cluCenter(numcluster,:) = mymean; %储存种子点 clusMemcell{numcluster} = mymembers; %存储类索引 cluvotes(:,numcluster) = thisvotes; %存储投票项 end break; end end initpind = find(~beenvisited); %更新点的序列，排序已经被循环的点 numn = length(initpind); %更新指针的长度，为0时结束循环 end [~,idx] = max(cluvotes,[],2); %根据投票项选择最多的聚类 if nargout > 2 clucell = cell(numcluster,1); for cn = 1:numcluster mymembers = find(idx == cn); clucell{cn} = mymembers; end end end