Matlab多光谱kmeans聚类分割

Matlab多光谱kmeans聚类分割

参考代码：https://blog.csdn.net/ma7856728/article/details/84891530
https://blog.csdn.net/qq_37970770/article/details/105606207

---

入驻csdn的小白 希望和大家一起学习进步 流水账似的有点凌乱 思绪还是稳得昂嘻嘻

---

上篇写到ENVI软件的分割工具出现错分的情况，按照之前所说的matlab设定阈值提取得到的图像没有地理意义，需要进行配准，并且只能得到每类的二值图，根据分类建立ROIs，提取值到点，总是提取到角点，不能得到像元值，如下图，而且arcgis和envi软件同时使用会遇到存储格式转换的问题，甚至在txt转excel时1600行就能拖卡电脑。
因此决定暂时不考虑envi、arcgis等软件，只用matlab进行处理。

笔者处理的影像为多光谱影像，matlab自带的kmeans函数可以处理多维数据，但每一维都要写成列向量的形式，通过查阅大家自己编写的函数以及不同的误差处理方式，了解到matlab自带的kmeans算法更优。

clc
clear
close all
I = imread('SPOT5多光谱.tif'); 
[M,N,L] = size(I);
%构造样本空间
A = reshape(I(:, :, 1), M*N, 1);    % 将RGB分量各转为kmeans使用的数据格式n行，一样一样本
B = reshape(I(:, :, 2), M*N, 1);
C = reshape(I(:, :, 3), M*N, 1);
D = reshape(I(:, :, 4), M*N, 1);
K = 4;
dat = [A B C D];  % 四个分量组成样本的特征，每个样本有4个属性值，共width*height个样本
c3 = kmeans(double(dat), K);    % 使用聚类算法分为K类
r3 = reshape(c3, M, N);     % 反向转化为图片形式
figure, imshow(label2rgb(r3))   % 显示分割结果
title('matlab库函数');

笔者的目的是分割每一类，因此编写了如下代码将结果导出

index = find(c3 == 1)%第i类 
temp1 = I(:,:,1)%第1个波段
temp2 = I(:,:,2)%第2个波段
temp3 = I(:,:,3)%第3个波段
temp4 = I(:,:,4)%第4个波段

b1= temp1(index)
b2= temp2(index)
b3= temp3(index)
b4= temp4(index)

class1=[b1,b2,b3,b4];
xlswrite('class1.xls',class1);   

为了方便后续研究的计算，需要修改代码，赋予地理坐标。

data=imread('D:\DDDDDDDDESK\异质性分析\样地100.tif');%读取纯数据

[multi_data,r]=geotiffread('D:\DDDDDDDDESK\异质性分析\样地100.tif');  % read the geo information

info=geotiffinfo('D:\DDDDDDDDESK\异质性分析\样地100.tif');   % read the geo information

[fl,s,b] = size(data);

dat = zeros(fl*s,b);%matlab K-means算法要求输入矩阵是一个列向量组成的矩阵，列数为波段数，每一列为fl*s的像元值

for i=1:b
    
  dat(:,i) = reshape(data(:,:,i),fl*s,1);

end

class_result = kmeans(dat,5);

out_data = reshape(class_result,fl,s);
figure, imshow(label2rgb(out_data))  % 显示聚类结果

geotiffwrite('geo_outimage.tif',out_data,r,'GeoKeyDirectoryTag',info.GeoTIFFTags.GeoKeyDirectoryTag);

index = find(class_result == 1)%第i类 
temp1 = data(:,:,1)%第1个波段
temp2 = data(:,:,2)%第2个波段
temp3 = data(:,:,3)%第3个波段
temp4 = data(:,:,4)%第4个波段

b1= temp1(index)
b2= temp2(index)
b3= temp3(index)
b4= temp4(index)

[m,n]=find(class_result == 1)  %寻找每类像元，得到m*1的矩阵
m1=fix(m/fl)+1; n1=rem(m,fl)+1;  %根据线性顺序转变为行列号
[a1,a2] = pix2map(info.RefMatrix,1,1)%左上角像素大地坐标
pixelscale=8; %输入单个像元大小
x=a1+pixelscale*(n1-1);%计算大地坐标
y=a2-pixelscale*(m1-1);
location=[x,y]; %double和uint8格式的数据存在不同的sheet中
class1=[b1,b2,b3,b4,b5];
xlswrite('100-class1.xls',class1,'Sheet1'); 
xlswrite('100-class1.xls',location,'Sheet2'); 

得到对应像元每个波段值和大地坐标，分别存储在同一表格的不同sheet下，随后进行每一类的抽样、变异函数分析、克里金插值等。
对变异函数分析后的样地进行比较range值与样地大小，某些情况需要分区进行kmeans处理，只需要改变初始左上角坐标和for循环中的像元个数即可。

得到的结果很明显有农田的规则格网状，统计分类结果如下：