MATLAB实现遥感图像分类——K均值算法

前言

这些代码均是使用最基础的方法，通过一步一步迭代过程来理解算法的原理及实现过程，并不采用于实用工程，读者以此作为学习参考即可。

非监督分类——k均值算法

图像分类的概念：将图像中的像素点划分为不同的类别，每个类别具有特定的物理含义。例如不同的地物类型，包括草地、水泥地、建筑、道路、森林、山地等。
图像分类的两个步骤：特征提取与分类算法。
特征提取：如何描述每一个像素点，颜色特征分量。
分类算法：K均值聚类算法。

监督分类与非监督分类

监督分类又称训练场地法、训练分类法，是以建立统计识别函数为理论基础、依据典型样本训练方法进行分类的技术，即根据已知训练区提供的样本，通过选择特征参数，求出特征参数作为决策规则，建立判别函数以对各待分类影像进行的图像分类。
无监督分类是以不同影像地物在特征空间中类别特征的差别为依据的一种无先验类别标准的图像分类，是以集群为理论基础，通过计算机对图像进行集聚统计分析的方法。根据待分类样本特征参数的统计特征，建立决策规则来进行分类。

K均值聚类算法

分类结果

由于采用rgb分量作为特征向量，分为三类，所以颜色接近的像素便分为一类，并在伪彩色处理中显示成相同的颜色。由于初始中心随机产生，所以表示三类的色彩也会随机变化，但分类结果均一致。

像素点分布的散点图

图中，蓝色的点为所有像素点的位置，红色的星号标注了第一类聚类中心的坐标位置变化，绿色的星号标注了第二类聚类中心的坐标位置变化，紫色的星号标注了第三类聚类中心的坐标位置变化

平均距离随迭代次数变化关系

这张图可以明显看出平均距离随着迭代次数逐渐缩小，在实验中存在减小后又突然增加的现象，推测是由于初始点的随机选择，导致存在某一类距离发生突变，而平均距离随着波动较大的类别距离发生了骤增，但从分类结果来看，各类别的像素点随着迭代次数明显集中。

源代码

clc;
clearvars
close all;

Im1 = imread('buildings91.tif');%读取图像
figure(1)
subplot(2,2,1),imshow(Im1);
ir=double(Im1(:,:,1));%分解图像红色分量
ig=double(Im1(:,:,2));%绿色分量
ib=double(Im1(:,:,3));%蓝色分量
[m,n]=size(ir);%取图像的大小
T = zeros(m,n);%设立类别标记表
d=zeros(3,1);%暂存到各类距离一维矩阵
dmin=zeros(256,1);%记录最短距离
cux=round(unifrnd(1,m,3,1));%随机产生三组坐标作为初始聚类中心
cuy=round(unifrnd(1,n,3,1));
Zz=zeros(n,9);%聚类中心保留矩阵
Zz(1,1)=ir(cux(1),cuy(1));
Zz(1,2)=ig(cux(1),cuy(1));
Zz(1,3)=ib(cux(1),cuy(1));
Zz(1,4)=ir(cux(2),cuy(2));
Zz(1,5)=ig(cux(2),cuy(2));
Zz(1,6)=ib(cux(2),cuy(2));
Zz(1,7)=ir(cux(3),cuy(3));
Zz(1,8)=ig(cux(3),cuy(3));
Zz(1,9)=ib(cux(3),cuy(3));%读取3个初始聚类中心的各色分量值
for count=1:256
    ddr1=0;ddg1=0;ddb1=0;%统计类别1中各分量的和
    ddr2=0;ddg2=0;ddb2=0;%统计类别2中各分量的和
    ddr3=0;ddg3=0;ddb3=0;%统计类别3中各分量的和
    d1=0;d2=0;d3=0;dmins=0;
    for i=1:m
        for j=1:n
            d(1)=sqrt((ir(i,j)-Zz(count,1))^2+(ig(i,j)-Zz(count,2))^2+(ib(i,j)-Zz(count,3))^2);%与第一类的距离
            d(2)=sqrt((ir(i,j)-Zz(count,4))^2+(ig(i,j)-Zz(count,5))^2+(ib(i,j)-Zz(count,6))^2);%与第二类的距离
            d(3)=sqrt((ir(i,j)-Zz(count,7))^2+(ig(i,j)-Zz(count,8))^2+(ib(i,j)-Zz(count,9))^2);%与第三类的距离
            dmins=dmins+min(d);
            if d(1)<d(2)&d(1)<d(3)%判断属于第一类
                T(i,j)=0;
                ddr1=ddr1+ir(i,j);
                ddg1=ddg1+ig(i,j);
                ddb1=ddb1+ib(i,j);
                d1=d1+1;
            elseif d(2)<d(1)&d(2)<d(3)%判断属于第二类
                T(i,j)=1;
                ddr2=ddr2+ir(i,j);
                ddg2=ddg2+ig(i,j);
                ddb2=ddb2+ib(i,j);
                d2=d2+1;
            else%判断属于第三类
                T(i,j)=2;
                ddr3=ddr3+ir(i,j);
                ddg3=ddg3+ig(i,j);
                ddb3=ddb3+ib(i,j);
                d3=d3+1;
            end
        end
    end
    dmin(count)=dmins/(m*n);
    Zz(count+1,1)=ddr1/d1;
    Zz(count+1,2)=ddg1/d1;
    Zz(count+1,3)=ddb1/d1;
    Zz(count+1,4)=ddr2/d2;
    Zz(count+1,5)=ddg2/d2;
    Zz(count+1,6)=ddb2/d2;
    Zz(count+1,7)=ddr3/d3;
    Zz(count+1,8)=ddg3/d3;
    Zz(count+1,9)=ddb3/d3;%计算各平均值
    if sum(Zz(count,:)==Zz(count+1,:))==9%判断是否结束循环
        break;
    end
end
figure(1)
subplot(2,2,2),imshow(T,[]);
tu1=zeros(m*n,3);%第一类点
tu2=zeros(m*n,3);%第二类点
tu3=zeros(m*n,3);%第三类点
%对灰度图进行伪彩色处理
[M N] = size(T);
T=T*255/(max(max(T)));
T2=zeros(M,N,3);                                 %初始化三通道
for x=1:M
   for y=1:N
      if T(x,y)<=127                                 % R
          T2(x,y,1)=0;
      elseif T(x,y)<=191
          T2(x,y,1)=4*T(x,y)-510;
      else
          T2(x,y,1)=255;
      end
      if T(x,y)<=63                                  % G
          T2(x,y,2)=254-4*T(x,y);
      elseif T(x,y)<=127
          T2(x,y,2)=4*T(x,y)-254;
      elseif T(x,y)<=191
          T2(x,y,2)=255;
      else
          T2(x,y,2)=1022-4*T(x,y);
      end
      if T(x,y)<=63                                 % B
          T2(x,y,3)=255;
      elseif T(x,y)<=127
          T2(x,y,3)=510-4*T(x,y);
      else
          T2(x,y,3)=0;
      end
   end
end
figure(1)
subplot(2,2,3),imshow(uint8(T2));
ir1=zeros(m*n,1);ig1=zeros(m*n,1);ib1=zeros(m*n,1);
k=1;
for i=1:m
    for j=1:n
        ir1(k)=ir(i,j);
        ig1(k)=ig(i,j);
        ib1(k)=ib(i,j);
        k=k+1;
    end
end
figure(2)
scatter3(ir1,ig1,ib1,'.');%生成散点图
xlabel('R');ylabel('G');zlabel('B');
hold on;scatter3(Zz(:,1),Zz(:,2),Zz(:,3),'*','r');
hold on;scatter3(Zz(:,4),Zz(:,5),Zz(:,6),'*','g');
hold on;scatter3(Zz(:,7),Zz(:,8),Zz(:,9),'*','m');
counts=1:count;
figure(3)
plot(counts,dmin(1:count));