《基于星座聚类的星敏感器导航星优选算法研究》的复现

简介

一直对公司采用的指向模型修正过程中的选星过程很感兴趣，想知道是否有相关的文献。

的确有找到相关文献，想就此进行复现。先对下文进行复现。

吴峰, 朱锡芳, 相入喜, 等. 基于星座聚类的星敏感器导航星优选算法研究[J]. 微电子学与计算机, 2018, 35(3): 140-144.

主要涉及Mean Shift（均值漂移）算法

---

起因

一般对原始星表恒星数据做星等阈值处理， 即可得到导航星星库． 为提高星图识别稳定性， 增强星敏感器系统性能， 通常要求导航星在天球上分布均匀． 为此， 国内外许多研究人员提出了进一步筛选导航星的算法．

---

相关文献：

---

算法过程

通常星座区域的恒星密度高， 星座和星座之间的区域星密度低． 为此， 本文首先借鉴Mean Shift（均值漂移）算法， 提出星座聚类算法， 将视场内的恒星划分到不同的星座． 然后， 提出基于星座聚类的导航星优选算法， 通过每次选择星数最多的一个星座， 并删除离该星座中心最近的一颗恒星， 最终实现导航星的均匀分布

Mean Shift（均值漂移）是基于密度的非参数聚类算法，其算法思想是假设不同簇类的数据集符合不同的概率密度分布，找到任一样本点密度增大的最快方向（最快方向的含义就是Mean Shift），样本密度高的区域对应于该分布的最大值，这些样本点最终会在局部密度最大值收敛，且收敛到相同局部最大值的点被认为是同一簇类的成员。（https://cloud.tencent.com/developer/article/1459530）

1.核密度估计

Mean Shift算法用核函数估计样本的密度，最常用的核函数是高斯核。它的工作原理是在数据集上的每一个样本点都设置一个核函数，然后对所有的核函数相加，得到数据集的核密度估计（kernel density estimation）。

假设我们有大小为n的d维数据集

如下图，我们用高斯核估计一维数据集的密度，每个样本点都设置了以该样本点为中心的高斯分布，累加所有的高斯分布，得到该数据集的密度。

其中虚线表示每个样本点的高斯核，实线表示累加所有样本高斯核后的数据集密度。因此，我们通过高斯核来得到数据集的密度。

5.图像分割示例

mean shift通过对像素空间进行聚类，达到图像分割的目的。

我们对下图进行图像分割：

Mean Shift算法的优缺点

优点：

不需要设置簇类的个数；

可以处理任意形状的簇类；

算法只需设置带宽这一个参数，带宽影响数据集的核密度估计

算法结果稳定，不需要进行类似K均值的样本初始化

缺点：

聚类结果取决于带宽的设置，带宽设置的太小，收敛太慢，簇类个数过多；带宽设置的太大，一些簇类可能会丢失。

对于较大的特征空间，计算量非常大。

图解过程

https://blog.csdn.net/google19890102/article/details/51030884

matlab代码

https://blog.csdn.net/HJ199404182515/article/details/121694298

https://www.cnblogs.com/kailugaji/p/11646167.html

函数代码为MeanShiftCluster.m

function [clustCent,data2cluster,cluster2dataCell] = MeanShiftCluster(dataPts,bandWidth,plotFlag)
%perform MeanShift Clustering of data using a flat kernel
%
% ---INPUT---
% dataPts           - input data, (numDim x numPts)
% bandWidth         - is bandwidth parameter (scalar)
% plotFlag          - display output if 2 or 3 D    (logical)
% ---OUTPUT---
% clustCent         - is locations of cluster centers (numDim x numClust)
% data2cluster      - for every data point which cluster it belongs to (numPts)
% cluster2dataCell  - for every cluster which points are in it (numClust)
%
% Bryan Feldman 02/24/06
% MeanShift first appears in
% K. Funkunaga and L.D. Hosteler, "The Estimation of the Gradient of a
% Density Function, with Applications in Pattern Recognition"
 
 
%*** Check input ****
if nargin < 2
    error('no bandwidth specified')
end
 
if nargin < 3
    plotFlag = true;
    plotFlag = false;
end
 
%**** Initialize stuff ***
[numDim,numPts] = size(dataPts);
numClust        = 0;
bandSq          = bandWidth^2;
initPtInds      = 1:numPts;
maxPos          = max(dataPts,[],2);                    %biggest size in each dimension
minPos          = min(dataPts,[],2);                    %smallest size in each dimension
boundBox        = maxPos-minPos;                        %bounding box size
sizeSpace       = norm(boundBox);                       %indicator of size of data space
stopThresh      = 1e-3*bandWidth;                       %when mean has converged
clustCent       = [];                                   %center of clust
beenVisitedFlag = zeros(1,numPts);                      %track if a points been seen already
numInitPts      = numPts;                               %number of points to posibaly use as initilization points
clusterVotes    = zeros(1,numPts);                      %used to resolve conflicts on cluster membership
 
 
while numInitPts
 
    tempInd         = ceil( (numInitPts-1e-6)*rand);        %pick a random seed point
    stInd           = initPtInds(tempInd);                  %use this point as start of mean
    myMean          = dataPts(:,stInd);                     % intilize mean to this points location
    myMembers       = [];                                   % points that will get added to this cluster                         
    thisClusterVotes    = zeros(1,numPts);                  %used to resolve conflicts on cluster membership
 
    while 1     %loop untill convergence
         
        sqDistToAll = sum((repmat(myMean,1,numPts) - dataPts).^2);    %dist squared from mean to all points still active
        inInds      = find(sqDistToAll < bandSq);                     %points within bandWidth
        thisClusterVotes(inInds) = thisClusterVotes(inInds)+1;        %add a vote for all the in points belonging to this cluster
         
         
        myOldMean   = myMean;                                   %save the old mean
        myMean      = mean(dataPts(:,inInds),2);                %compute the new mean
        myMembers   = [myMembers inInds];                       %add any point within bandWidth to the cluster
        beenVisitedFlag(myMembers) = 1;                         %mark that these points have been visited
         
        %*** plot stuff ****
        if plotFlag
            figure(1),clf,hold on
            if numDim == 2
                plot(dataPts(1,:),dataPts(2,:),'.')
                plot(dataPts(1,myMembers),dataPts(2,myMembers),'ys')
                plot(myMean(1),myMean(2),'go')
                plot(myOldMean(1),myOldMean(2),'rd')
                pause
            end
        end
 
        %**** if mean doesn't move much stop this cluster ***
        if norm(myMean-myOldMean) < stopThresh
             
            %check for merge posibilities
            mergeWith = 0;
            for cN = 1:numClust
                distToOther = norm(myMean-clustCent(:,cN));     %distance from posible new clust max to old clust max
                if distToOther < bandWidth/2                    %if its within bandwidth/2 merge new and old
                    mergeWith = cN;
                    break;
                end
            end
             
             
            if mergeWith > 0    % something to merge
                clustCent(:,mergeWith)       = 0.5*(myMean+clustCent(:,mergeWith));              %record the max as the mean of the two merged (I know biased twoards new ones)
                %clustMembsCell{mergeWith}    = unique([clustMembsCell{mergeWith} myMembers]);   %record which points inside
                clusterVotes(mergeWith,:)    = clusterVotes(mergeWith,:) + thisClusterVotes;     %add these votes to the merged cluster
            else    %its a new cluster
                numClust                    = numClust+1;                    %increment clusters
                clustCent(:,numClust)       = myMean;                        %record the mean 
                %clustMembsCell{numClust}    = myMembers;                    %store my members
                clusterVotes(numClust,:)    = thisClusterVotes;
            end
 
            break;
        end
 
    end
     
     
    initPtInds      = find(beenVisitedFlag == 0);           %we can initialize with any of the points not yet visited
    numInitPts      = length(initPtInds);                   %number of active points in set
 
end
 
[val,data2cluster] = max(clusterVotes,[],1);                %a point belongs to the cluster with the most votes
 
%*** If they want the cluster2data cell find it for them
if nargout > 2
    cluster2dataCell = cell(numClust,1);
    for cN = 1:numClust
        myMembers = find(data2cluster == cN);
        cluster2dataCell{cN} = myMembers;
    end
end

9.32   8.16  
8.45    -4.4
7.67    0.73
7.3    -3.38
5.72  2.38
5.5  -3.7
5.31  -9.36
4.31  -5.68
4.08  -6.6
2.52  3.26
1.24 -3.63
0.24 7.19
-0.45  4.09
-0.76  -2.59
-3.54  8.17
-4.49  0.34
-4.61  9.58
-5.76  -6.98
-6.55  7.11
-7.79  8.99
-9.40  -3.37

测试代码testMeanShift.m为

clear
clc
profile on
 
bandwidth = 2;
%% 加载数据
data_load=dlmread('Copy_of_gauss_data.txt');
[~,dim]=size(data_load);
data=data_load;
x=data';
%% 聚类
tic
[clustCent,point2cluster,clustMembsCell] = MeanShiftCluster(x,bandwidth);
% clustCent：聚类中心 D*K, point2cluster：聚类结果 类标签, 1*N
toc
%% 作图
numClust = length(clustMembsCell);
figure(2),clf,hold on
cVec = 'bgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmyk';%, cVec = [cVec cVec];
for k = 1:min(numClust,length(cVec))
    myMembers = clustMembsCell{k};
    myClustCen = clustCent(:,k);
%     plot(x(1,myMembers),x(2,myMembers),[cVec(k) '.']);
    plot(x(2,myMembers),x(1,myMembers),[cVec(k) '.']);
    plot(myClustCen(2),myClustCen(1),'o','MarkerEdgeColor','k','MarkerFaceColor',cVec(k), 'MarkerSize',10)

%     plot(myClustCen(1),myClustCen(2),'o','MarkerEdgeColor','k','MarkerFaceColor',cVec(k), 'MarkerSize',10)
end
for j=1:size(data,1)
    text(data(j,2)+0.2,data(j,1)-0.2,num2str(j));
end

title(['no shifting, numClust:' int2str(numClust)])

结果不唯一，有时为18类。

clustCent 2*18 为中心坐标

point2cluster 1*21 为类别

clustMembsCell =

18×1 cell 数组

{[ 19]}
{[ 11]}
{[15 17]}
{[ 3]}
{[ 21]}
{[ 10]}
{[ 13]}
{[ 1]}
{[ 16]}
{[ 14]}
{[2 4 6]}
{[ 8 9]}
{[ 7]}
{[ 18]}
{[ 12]}
{[ 5]}
{[ 20]}

根据文章算法修改后MeanShiftCluster.m

function [clustCent,data2cluster,cluster2dataCell] = MeanShiftCluster(dataPts,bandWidth,plotFlag)
%perform MeanShift Clustering of data using a flat kernel
%
% ---INPUT---
% dataPts           - input data, (numDim x numPts)
% bandWidth         - is bandwidth parameter (scalar)
% plotFlag          - display output if 2 or 3 D    (logical)
% ---OUTPUT---
% clustCent         - is locations of cluster centers (numDim x numClust)
% data2cluster      - for every data point which cluster it belongs to (numPts)
% cluster2dataCell  - for every cluster which points are in it (numClust)
%
% Bryan Feldman 02/24/06
% MeanShift first appears in
% K. Funkunaga and L.D. Hosteler, "The Estimation of the Gradient of a
% Density Function, with Applications in Pattern Recognition"
 
 
%*** Check input ****
if nargin < 2
    error('no bandwidth specified')
end
 
if nargin < 3
    plotFlag = true;
    plotFlag = false;
end

%**** Initialize stuff ***
[numDim,numPts] = size(dataPts);
numClust        = 0;
bandSq          = bandWidth^2;
initPtInds      = 1:numPts;
maxPos          = max(dataPts,[],2);                    %biggest size in each dimension
minPos          = min(dataPts,[],2);                    %smallest size in each dimension
boundBox        = maxPos-minPos;                        %bounding box size
sizeSpace       = norm(boundBox);                       %indicator of size of data space
stopThresh      = 1e-3*bandWidth;                       %when mean has converged
clustCent       = [];                                   %center of clust
beenVisitedFlag = zeros(1,numPts);                      %track if a points been seen already
numInitPts      = numPts;                               %number of points to posibaly use as initilization points
clusterVotes    = zeros(1,numPts);                      %used to resolve conflicts on cluster membership
nowNum = 1;
outCom = zeros(1,numPts);
 
while numInitPts

    tempInd         = ceil( (numInitPts-1e-6)*rand);        %pick a random seed point
%     stInd           = initPtInds(tempInd);                  %use this point as start of mean
    stInd           = min(initPtInds);                  %use this point as start of mean
    myMean          = dataPts(:,stInd);                     % intilize mean to this points location
    myMembers       = [];                                   % points that will get added to this cluster                         
    thisClusterVotes    = zeros(1,numPts);                  %used to resolve conflicts on cluster membership
%     beyondMe=[];
    while 1     %loop untill convergence
         
        sqDistToAll = sum((repmat(myMean,1,numPts) - dataPts).^2);    %dist squared from mean to all points still active
        inInds      = find(sqDistToAll < bandSq);                     %points within bandWidth

        if (length(inInds)>1 && min(outCom(inInds))==0)% 2，6
            beyondMe = inInds;
            beyondMe(beyondMe == initPtInds(1)) = [];
            if(find(outCom(beyondMe) >0))
                thisCluster = min(outCom(find(outCom(beyondMe) >0)));
                outCom(stInd)=outCom(beyondMe(find(outCom(beyondMe) >0)));
                
            else % 4
                outCom(stInd)=nowNum;
                outCom(beyondMe)=nowNum;
                nowNum = nowNum + 1;
            end
        elseif length(inInds)==1 % 1、3、5
            outCom(stInd)=nowNum;
            nowNum = nowNum + 1;
        end
        thisClusterVotes(inInds) = thisClusterVotes(inInds)+1;        %add a vote for all the in points belonging to this cluster
        myOldMean   = myMean;                                   %save the old mean
        myMean      = mean(dataPts(:,inInds),2);                %compute the new mean
        myMembers   = [myMembers inInds];                       %add any point within bandWidth to the cluster
        beenVisitedFlag(myMembers) = 1;                         %mark that these points have been visited
 
        %**** if mean doesn't move much stop this cluster ***
        if norm(myMean-myOldMean) < stopThresh
             
            %check for merge posibilities
            mergeWith = 0;
            for cN = 1:numClust
                distToOther = norm(myMean-clustCent(:,cN));     %distance from posible new clust max to old clust max
                if distToOther < bandWidth/2                    %if its within bandwidth/2 merge new and old
                    mergeWith = cN;
                    break;
                end
            end
             
             
            if mergeWith > 0    % something to merge
                clustCent(:,mergeWith)       = 0.5*(myMean+clustCent(:,mergeWith));              %record the max as the mean of the two merged (I know biased twoards new ones)
                %clustMembsCell{mergeWith}    = unique([clustMembsCell{mergeWith} myMembers]);   %record which points inside
                clusterVotes(mergeWith,:)    = clusterVotes(mergeWith,:) + thisClusterVotes;     %add these votes to the merged cluster
            else    %its a new cluster
                numClust                    = outCom(stInd);                    %increment clusters
                clustCent(:,numClust)       = myMean;                        %record the mean 
                %clustMembsCell{numClust}    = myMembers;                    %store my members
                clusterVotes(numClust,:)    = thisClusterVotes;
            end

            break;
        end
 
    end
     
    
    initPtInds      = find(beenVisitedFlag == 0);           %we can initialize with any of the points not yet visited
    numInitPts      = length(initPtInds);                   %number of active points in set
    [val,data2cluster] = max(clusterVotes,[],1);
end

data2cluster = outCom;

%*** If they want the cluster2data cell find it for them
if nargout > 2
    cluster2dataCell = cell(numClust,1);
    for cN = 1:numClust
        myMembers = find(data2cluster == cN);
        cluster2dataCell{cN} = myMembers;
    end
end

outCom =

 1     2     3     2     4     2     5     6     6     7     8     9    10    11    12    13    12    14    15    16    17

clustCent 2*17 为中心坐标

point2cluster 1*21 为类别 1 1 3 2 4 2 5 6 6 7 8 9 10 11 12 13 12 14 15 16 17

根据文章算法修改后MeanShiftCluster.m

clustMembsCell =

18×1 cell 数组

{[ 19]}
{[ 11]}
{[15 17]}
{[ 3]}
{[ 21]}
{[ 10]}
{[ 13]}
{[ 1]}
{[ 16]}
{[ 14]}
{[2 4 6]}
{[ 8 9]}
{[ 7]}
{[ 18]}
{[ 12]}
{[ 5]}
{[ 20]}

经过第一轮循环，删除了第4/8/15个星。

再正式修改为MeanShiftCluster.m

function [clustCent,data2cluster,cluster2dataCell] = MeanShiftCluster(dataPts,bandWidth,plotFlag)
%perform MeanShift Clustering of data using a flat kernel
%
% ---INPUT---
% dataPts           - input data, (numDim x numPts)
% bandWidth         - is bandwidth parameter (scalar)
% plotFlag          - display output if 2 or 3 D    (logical)
% ---OUTPUT---
% clustCent         - is locations of cluster centers (numDim x numClust)
% data2cluster      - for every data point which cluster it belongs to (numPts)
% cluster2dataCell  - for every cluster which points are in it (numClust)
%
% Bryan Feldman 02/24/06
% MeanShift first appears in
% K. Funkunaga and L.D. Hosteler, "The Estimation of the Gradient of a
% Density Function, with Applications in Pattern Recognition"
 
 
%*** Check input ****
if nargin < 2
    error('no bandwidth specified')
end
 
if nargin < 3
    plotFlag = true;
    plotFlag = false;
end

%**** Initialize stuff ***
[numDim,numPts] = size(dataPts);
numClust        = 0;
bandSq          = bandWidth^2;
initPtInds      = 1:numPts;
maxPos          = max(dataPts,[],2);                    %biggest size in each dimension
minPos          = min(dataPts,[],2);                    %smallest size in each dimension
boundBox        = maxPos-minPos;                        %bounding box size
sizeSpace       = norm(boundBox);                       %indicator of size of data space
stopThresh      = 1e-3*bandWidth;                       %when mean has converged
clustCent       = [];                                   %center of clust
beenVisitedFlag = zeros(1,numPts);                      %track if a points been seen already
numInitPts      = numPts;                               %number of points to posibaly use as initilization points
clusterVotes    = zeros(1,numPts);                      %used to resolve conflicts on cluster membership
nowNum = 1;
outCom = zeros(1,numPts);
 
while numInitPts

    tempInd         = ceil( (numInitPts-1e-6)*rand);        %pick a random seed point
    stInd           = initPtInds(tempInd);                  %use this point as start of mean
    myMean          = dataPts(:,stInd);                     % intilize mean to this points location
    myMembers       = [];                                   % points that will get added to this cluster                         
    thisClusterVotes    = zeros(1,numPts);                  %used to resolve conflicts on cluster membership
%     beyondMe=[];
    while 1     %loop untill convergence
         
        sqDistToAll = sum((repmat(myMean,1,numPts) - dataPts).^2);    %dist squared from mean to all points still active
        inInds      = find(sqDistToAll < bandSq);                     %points within bandWidth

        if (length(inInds)>1 && min(outCom(inInds))==0)% 2，6
            beyondMe = inInds;
            beyondMe(beyondMe == stInd) = [];
            if(find(outCom(beyondMe) >0))
                thisCluster = min(outCom(find(outCom(beyondMe) >0)));
                outCom(stInd)=min(outCom(beyondMe(find(outCom(beyondMe) >0))));
                outCom(beyondMe)=min(outCom(beyondMe(find(outCom(beyondMe) >0))));
            else % 4
                outCom(stInd)=nowNum;
                outCom(beyondMe)=nowNum;
                nowNum = nowNum + 1;
            end
        elseif length(inInds)==1 % 1、3、5
            outCom(stInd)=nowNum;
            nowNum = nowNum + 1;
        end
        thisClusterVotes(inInds) = thisClusterVotes(inInds)+1;        %add a vote for all the in points belonging to this cluster
        myOldMean   = myMean;                                   %save the old mean
        myMean      = mean(dataPts(:,inInds),2);                %compute the new mean
        myMembers   = [myMembers inInds];                       %add any point within bandWidth to the cluster
        beenVisitedFlag(myMembers) = 1;                         %mark that these points have been visited
 
        %**** if mean doesn't move much stop this cluster ***
        if norm(myMean-myOldMean) < stopThresh
             
            %check for merge posibilities
            mergeWith = 0;
            for cN = 1:numClust
                distToOther = norm(myMean-clustCent(:,cN));     %distance from posible new clust max to old clust max
                if distToOther < bandWidth/2                    %if its within bandwidth/2 merge new and old
                    mergeWith = cN;
                    break;
                end
            end
             
             
            if mergeWith > 0    % something to merge
                clustCent(:,mergeWith)       = 0.5*(myMean+clustCent(:,mergeWith));              %record the max as the mean of the two merged (I know biased twoards new ones)
                %clustMembsCell{mergeWith}    = unique([clustMembsCell{mergeWith} myMembers]);   %record which points inside
                clusterVotes(mergeWith,:)    = clusterVotes(mergeWith,:) + thisClusterVotes;     %add these votes to the merged cluster
            else    %its a new cluster
                numClust                    = outCom(stInd);                    %increment clusters
                clustCent(:,numClust)       = myMean;                        %record the mean 
                %clustMembsCell{numClust}    = myMembers;                    %store my members
                clusterVotes(numClust,:)    = thisClusterVotes;
            end

            break;
        end
 
    end
     
    
    initPtInds      = find(beenVisitedFlag == 0);           %we can initialize with any of the points not yet visited
    numInitPts      = length(initPtInds);                   %number of active points in set
    [val,data2cluster] = max(clusterVotes,[],1);
end

data2cluster = outCom;

%*** If they want the cluster2data cell find it for them
if nargout > 2
    cluster2dataCell = cell(numClust,1);
    for cN = 1:numClust
        myMembers = find(data2cluster == cN);
        cluster2dataCell{cN} = myMembers;
    end
end

testMeanShift.m为

clear
clc
profile on
 
bandwidth = 2;
%% 加载数据
data_load=dlmread('Copy_of_gauss_data.txt');
[~,dim]=size(data_load);
data=data_load;
x=data';
%% 聚类
tic

    while length(data)>6
        [clustCent,point2cluster,clustMembsCell] = MeanShiftCluster(data.',bandwidth);
        numClust = length(clustMembsCell);
        while (length(find(point2cluster==mode(point2cluster)))>1)
            [val,idx]=min(dist(clustCent(:,mode(point2cluster)).',data.'))   ;
            data(idx,:)=[];
            x(:,idx)=[];
            [clustCent,point2cluster,clustMembsCell] = MeanShiftCluster(data.',bandwidth);
        end
        bandwidth = bandwidth + 2;
    end

figure(2),clf,hold on;axis equal
cVec = 'bgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmyk';%, cVec = [cVec cVec];
for k = 1:min(numClust,length(cVec))
    myMembers = clustMembsCell{k};
    myClustCen = clustCent(:,k);
    plot(x(2,myMembers),x(1,myMembers),[cVec(k) '.']);
    plot(myClustCen(2),myClustCen(1),'o','MarkerEdgeColor','k','MarkerFaceColor',cVec(k), 'MarkerSize',10)
end
title(['no shifting, numClust:' int2str(numClust)])
hold on 
toc

得到不多于6个点的几个测试结果为：

通常以玻尔兹曼熵评价导航星的全天球分布均匀性， 以视场内捕获恒星数的最大值、 最小值、 均值和标准差等评价局部天球分布均匀。

但我还不知道如何求解玻尔兹曼熵。（下文求解）

---

带入星图数据（下代码有误，请看后一篇博客）

clear
clc
profile on
load("附件2  简易星表.mat");
bandwidth = 5;
%% 加载数据
data_load=dlmread('Copy_of_gauss_data.txt');
[~,dim]=size(data_load);
data = star_data(:,2:3);
% data = data_load(:,1:2);
x=data';
%% 聚类
tic
    while length(data)>150
        [clustCent,point2cluster,clustMembsCell] = MeanShiftCluster(data.',bandwidth);
        numClust = length(clustMembsCell);
        while (length(find(point2cluster==mode(point2cluster)))>1)
            [val,idx]=min(dist(clustCent(:,mode(point2cluster)).',data.'))   ;
            data(idx,:)=[];
            x(:,idx)=[];
            [clustCent,point2cluster,clustMembsCell] = MeanShiftCluster(data.',bandwidth);
        end
        numClust = length(clustMembsCell);
        bandwidth = bandwidth + 2;
    end

figure(2),clf,hold on;axis equal
cVec = 'bgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmykbgrcmyk';%, cVec = [cVec cVec];
for k = 1:min(numClust,length(cVec))
    myMembers = clustMembsCell{k};
    myClustCen = clustCent(:,k);
    plot(x(2,myMembers),x(1,myMembers),[cVec(k) '.']);
    plot(myClustCen(2),myClustCen(1),'o','MarkerEdgeColor','k','MarkerFaceColor',cVec(k), 'MarkerSize',10)
end

title(['no shifting, numClust:' int2str(numClust)])
hold on 
toc

历时 3101.811945 秒。