中国地质大学武汉模式识别第二次试验kmeans

1、模式识别第二次实验 学号： 20121001873 班学号：07512213 姓名： 吴泽光 指导老师：马丽 中国地质大学（武汉）机械与电子信息学院 2015年4月题目：基于k-means算法的遥感图像分类实验目标：1. 掌握K-means算法原理2. 用MATLAB实现k-means算法，并进行结果分析实验内容：1. 采用zy3sample数据（400*400*4），该数据包含4个类别：土地；农作物/植被；大棚；道路。2. 关于初始聚类中心的选择，如果是采用随机选择方法，那么至少随机选择3次，并进行结果比较。（注意：道路占有的像素很少，所以随机选择4个聚类中心的话，很可能选择不到这种地物。可

2、以思考该如何解决这个问题？可以简单的先从带有标签的CMap（也就是自己做的ROI保存的图像）中选择，也就是当选择的4个聚类中心都正确的分别属于一类地物的话，看看k-means算法的结果，是否比在整个图像中随机选择要好。）3. 关于迭代终止条件。可以参考下面的一种方法：（1）设置一个阈值，当前后两次聚类中心的差别小于这个阈值的时候，停止迭代。（2）设置一个阈值，当前后两次聚类中所包含的点的数目变化小于这个阈值的时候，停止迭代。（3）设置最多迭代次数，也就是不管迭代有没有收敛，只要达到这个迭代次数就停止。4. 关于距离度量方法，可以选用欧氏距离，光谱角距离（提供代码），如果两个都选择，那么比较一下

3、聚类效果哪个好，或者也可以查资料采用其他更好的方法。5. 画出整个图像的聚类结果，每种类别采用一种彩色ujuj（提供ColorTable和画图代码）。基本要求：1.每位同学必须掌握k-mean算法的原理，以及会用MATLAB把算法实现，对图像进行分类。发挥部分：1. 关于初始聚类中心的选择，进行方法的比较和思考。2. 关于距离度量方法，进行方法的比较和思考。3. 实验的目的就是得到更好的聚类效果，不管采用什么方法，只要达到这个目的就可以。试验过程：1，k-means算法原理： K-means算法是很典型的基于距离的聚类算法，采用距离作为相似性的评价指标，即认为两个对象的距离越近，其相似度就越大

4、。该算法认为簇是由距离靠近的对象组成的，因此把得到紧凑且独立的簇作为最终目标。 算法过程如下：1)从N个数据中随机选取K个点作为类心2）对剩余的每个数据点测量其到每个类心的距离，并把它归到最近的类心的类3）重新计算已经得到的各个类的类心4）迭代23步直至新的类心与原类心相等或小于指定阈值，算法结束2，算法代码：1)首先选择每个类的平均值作为每个类的类心，类心迭代5次：% 读遥感数据，以及读入带标签的数据%在CMAP中求平均值，欧式距离,定阈值clear % read the data zy3sample1Img = multibandread('zy3sample1',400,

5、400,4,'float',0,'bsq','n','Band','Direct',1:4);%将高光谱数据读入GT = multibandread('zys',400,400,1,'uint8',0,'bsq','n','Band','Direct',1:1);%将ground truth读入（也就是带标签的CMap或说ROI信息读入）test_class=1:4;C=length(test_class);NbRow=

6、400;NbCol=400;NbDim=4;dataname='zy3'X=;Y=;A=zeros(C,4);for i=1:C row,col=find(GT=i);%总共21837个点有值 for j=1:length(row) data=reshape(Img(row(j),col(j),:),1,NbDim); A(i,:)=A(i,:)+data; X=X;data;%采样4列 Y=Y;i;%标签信息 end A(i,:)=A(i,:)./length(row);%平均聚类中心enddata=reshape(Img,160000,NbDim);%将原始数据变成2维数据

7、j=1;for i=1:C %求各点到类心的距离 D=Euclidian_distance(data,A(i,:); D1(:,j)=D; j=j+1;endc,d=min(D1,2);%返回每行的最小值的坐标v=;A1=zeros(4,4);D=;K-means部分：for n=1:5 A2=A1; for i=1:C %重新分类 N=; l=find(d=i); %寻找类别的地址 N=data(l,:); A1(i,:)=mean(N);%新的平均聚类中心,求每一列的均值 end for i=1:C %求各点到类心的距离 D=Euclidian_distance(data,A1(i,:);

8、%求出距离再次分类 D2(:,i)=D; end c,d=min(D2,2);%返回每行的最小值的坐标 end DATA=d;LINES=400;SAMPLES=400;ColorTable=160,82,45;% c1:soil 255, 255, 255 ;%c3:house 0 255 0;%c4 road 255,0,0 ;%c2:grassCOLORTABLE=ColorTable; ClassificationMap = GenerateClassificationMap( LINES, SAMPLES, COLORTABLE, DATA );imshow(Classificati

9、onMap);title('迭代5次后的分类结果')2）随机选择类心，对比程序：A=zeros(4,4);for i=1:4r=randperm(16000);r=r(1);A(i,:)=data(r,:);%随机数据做为类心end3）求数据之间的光谱角距离：%光谱角距离function D=SAM(data1,data2) ;% input：nsamp*nband%data1，data2,(n*d)的列向量为data，行向量为dimensions%data1未知光谱，data2已知光谱%结果：0到pi/2，越相似，D越小，越接近0.if nargin<2, %如果只有1

10、个输入，就是求data1之间的SAM距离 length_1=sqrt(sum(data1.*data1)'); % length_1：1*N; data:N*L D=real(acos(data1*data1'./(length_1'*length_1); D = min(D,D'); %保证对称 for i=1:size(D,1) D(i,i)=0; end %保证对角线值为0else %如果2个输入，就是求data1和data2之间的SAM距离 length_1=sqrt(sum(data1.*data1)'); % 1*N; data:N*L le

11、ngth_2=sqrt(sum(data2.*data2)'); D=real(acos(data1*data2'./(length_1'*length_2);endfunction D=Euclidian_distance(data1,data2) %data:N*Lif nargin<2 %求data1数据之间的欧氏距离阵D length_1=sum(data1.*data1)'); % 1*N; data:N*L ab=data1*data1' D=repmat(length_1,size(length_1,2),1) + repmat(le

12、ngth_1',1 size(length_1,2) -2*ab; D = min(D,D'); D=sqrt(D); for i=1:size(D,1) D(i,i)=0; endelse %求data1和data2数据之间的欧氏距离阵D length_1=sum(data1.*data1)',1); % 1*N; data:N*L length_2=sum(data2.*data2)',1); ab=data1*data2' data1=;data2=; D=repmat(length_1',1,size(length_2,2) + repmat(length_2,size(length_1,2),1) -2*ab; D=sqrt(D);endD=real(D);4）求数据之间的欧氏距离：3）算法结果：欧氏距离，以每个类的平均值为类心迭代5次后的结果：对比上个算法，随机选择类心的结果： 对比第一个算法，随机选择类心迭代一次后的结果： 对比第一个算法，选择平均类心循环一个后的结果： 对比上一个算法，选择光谱角后的结果： 4，结果分析： 总体的分类效果都不是很好。从迭代次数来看，选择每一类的平均值作为类心且只迭代一次的效果最好。随机选择类心分