效果图
Kmeans主函数
function Kmeans(data,kinds)
%传入一个图像矩阵进行非监督分类,根据光谱特征进行分类
%其中kinds为待分类的类别
[row,col,bands]=size(data);%求得多光谱的尺寸
%接下来为所选的类别产生随机数，但是由于只是熟悉算法，暂时不考虑精确性
%先进行去极化线性拉伸之后再分类
data1=zeros(size(data));
for i=1:bands
data1(:,:,i)=LineExtension(data(:,:,i));
end
%随机产生kinds个随机点作为聚类中心
clustercenter1=RandomPoints(kinds,row,col);%均值聚类向量中心点
temp=size(clustercenter1);%求临时尺寸
clustercenter=zeros(temp(1),bands);%聚类中心的平均值
for i=1:temp(1)
center_row=clustercenter1(i,1);%获得第i个聚类中心的行
center_col=clustercenter1(i,2);%获得第i个聚类中心的列
clustercenter(i,:)=data(center_row,center_col,:);%聚类中心像素的值应该是3行7列之类的
end
clustercenter_pre=zeros(size(clustercenter));%之前的聚类中心点
%由于最后展示的是分类的效果，是一种单波段的二维平面图
label=zeros(row,col);
threshold=0.01;
%选取好聚类中心之后就进行循环,选定聚类标准进行聚类---这里使用最短距离的方法进行聚类分析
while VectorDis(clustercenter,clustercenter_pre)>threshold
%当两次聚类中心的差距超过0.1就继续进行循环分类
clustercenter_pre=clustercenter;%将前一步的迭代更新
[clustercenter,label]=MinDis(data,clustercenter);%根据最短距离的方法进行非监督分类
end
%将计算出来的label标签使用imagesec进行显示
figure
imagesc(label);title('classified image');
end
%%
%通过欧式距离的方法进行重分类
function [clustercenter,label]=MinDis(data,clustercenter)
%声明一个标签数组用于存放分类的情况，clustercenter为聚类中心的平均值
[row,col,bands]=size(data);
label=zeros(row,col);
[kinds,~]=size(clustercenter);%只接受第一个参数，知道用户想要把图像分为多少类
%声明一个数组用于存放放进每个类别的总和，这样子就可以不用在每次收录进来的时候遍历图像数组
all_sum=clustercenter;
%声明一个数组用来统计各个聚类中心的数目
count_sum=ones(kinds,1);
%第一步遍历图像中的每一个像素位置的空间点
for i=1:row
for j=1:col
%遍历二维图像矩阵并计算该像素点与随机生成的各个聚类中心的距离
tempdata=reshape(data(i,j,:),[1,bands]);
mindis=VectorDis(double(tempdata),double(clustercenter(1,:)));
pixel_label=1;%先假设这个像素点距离第一类聚类中心最近
for k=2:kinds
%先将data(i,j,:)reshape
tempdis=VectorDis(double(tempdata),double(clustercenter(k,:)));
if mindis>tempdis
mindis=tempdis;
pixel_label=k;%表示将像素归于第k类
end
end
label(i,j)=pixel_label;%将标签数组归为第k类
%每收录一个点就重新计算聚类中心
%为了增加效率，只重新计算第k类的中心
all_sum(pixel_label,:)=tempdata + all_sum(pixel_label,:);
%将各个聚类中心的数目加一
count_sum(pixel_label,1)=count_sum(pixel_label,1)+1;
%重新将第pixel_label类计算聚类中心
clustercenter(pixel_label,:)= double(all_sum(pixel_label,:))/count_sum(pixel_label,1);
end
end
%根据归类后的状态，重新计算聚类中心
for i=1:kinds
sum=zeros(1,bands);
for j=1:row
for k=1:col
if label(j,k)==i
%如果当前像素点归为第i类
temp_=reshape(data(j,k,:),[1,bands]);
sum=sum+double(temp_);
end
end
end
%求得上述值之后就进行均值的重新赋值
clustercenter(i,:)=sum/(row*col);
end
end
%%
%这个函数算得二维范数的一个距离
function [dis]=VectorDis(vector1,vector2)
[row1,col1]=size(vector1);
sum1=0;
try
%判断用户传入的是行向量还是列向量
mainD=1;
if row1<col1
mainD=2;
end
if mainD==1
temp=vector2-vector1;
for i=1:row1
sum1=temp(i,1)^2+sum1;
end
%为了提高运算速率就不进行根号求解
dis=sqrt(sum1);%还是进行根号的求解
else
temp=vector2-vector1;
for i=1:col1
sum1=temp(1,i)^2+sum1;
end
%为了提高运算速率就不进行根号求解
dis=sqrt(sum1);%还是进行根号的求解
end
catch
error('the size of two vector you input must be the same size!');
end
end
%%
%下面的函数是生成随机坐标
function [randompoints]=RandomPoints(kinds,row,col)
randompoints=zeros(kinds,2);%就这个例子来说直接就是生成7行2列的数组进行坐标值的存储
for i=1:kinds
while true
flag=false;
%随机生成行坐标
m=randi(row,1);
%随机生成列坐标
n=randi(col,1);
%遍历已经有的坐标看是否已经存在，如果已经存在就进行相关的重生成操作
for j=1:kinds
if m==randompoints(j,1)&&n==randompoints(j,2)
flag=true;
break;
end
end
if flag is true
%do nothing
else
randompoints(i,1)=m;
randompoints(i,2)=n;
break;%跳出while循环
end
end
end
end
读取图像的函数
function data=read_ENVIimage3(imagefile)
feature('DefaultCharacterSet', 'UTF8');
%读取.img格式的二五年间，前提是图像以',img'后缀名
if length(imagefile)>=4 %如果位于同级目录下的长度大于4才进行操作
%主要的思路就是利用提高的img格式的文件寻找位于同级目录下的hdr文件
switch strcmp(imagefile(length(imagefile)-3:end),'dat')
case 0
hdrfilename=strcat(imagefile(1:(length(imagefile)-4)),'.hdr');
end
else
hdrfilename=strcat(imagefile,'hdr');
end
% end
% 以只读的方式打开头文件
fid = fopen(hdrfilename,'r');
info=fread(fid,'char=>char');
info=info';%转化为行向量
fclose(fid);
%查找列数
a=strfind(info,'samples = ');
b=length('samples = ');
c=strfind(info,'lines');
samples=[];
for i=a+b:c-1
samples=[samples,info(i)];
end
samples=str2num(samples);
%查找行数
a=strfind(info,'lines = ');
b=length('lines = ');
c=strfind(info,'bands');
lines=[];
for i=a+b:c-1
lines=[lines,info(i)];
end
lines=str2num(lines);
%查找波段数
a=strfind(info,'bands = ');
b=length('bands = ');
c=strfind(info,'data type');
bands=[];
for i=a+b:c-1
bands=[bands,info(i)];
end
bands=str2num(bands);
%查找数据类型
a=strfind(info,'data type = ');
b=length('data type = ');
c=strfind(info,'interleave');
datatype=[];
for i=a+b:c-1
datatype=[datatype,info(i)];
end
datatype=str2num(datatype);
precision=[];
switch datatype
case 1
precision='uint8=>uint8';%头文件中datatype=1对应ENVI中数据类型为Byte，对应MATLAB中数据类型为uint8
case 2
precision='int16=>int16';%头文件中datatype=2对应ENVI中数据类型为Integer，对应MATLAB中数据类型为int16
case 12
precision='uint16=>uint16';%头文件中datatype=12对应ENVI中数据类型为Unsighed Int，对应MATLAB中数据类型为uint16
case 3
precision='int32=>int32';%头文件中datatype=3对应ENVI中数据类型为Long Integer，对应MATLAB中数据类型为int32
case 13
precision='uint32=>uint32';%头文件中datatype=13对应ENVI中数据类型为Unsighed Long，对应MATLAB中数据类型为uint32
case 4
precision='float32=>float32';%头文件中datatype=4对应ENVI中数据类型为Floating Point，对应MATLAB中数据类型为float32
case 5
precision='double=>double';%头文件中datatype=5对应ENVI中数据类型为Double Precision，对应MATLAB中数据类型为double
otherwise
error('invalid datatype');%除以上几种常见数据类型之外的数据类型视为无效的数据类型
end
%查找数据格式
a=strfind(info,'header offset =');
c=strfind(info,'byte order ');
b=length('header offset = ');
headeroffset=[];
for i=a+b:c-1
headeroffset=[headeroffset,info(i)];
end
if headeroffset~='0'
headeroffset=strtrim(headeroffset);%删除字符串中前后的空格
end
%查找数据格式
a=strfind(info,'interleave = ');
b=length('interleave = ');
c=strfind(info,'file type');
interleave=[];
for i=a+b:c-1
interleave=[interleave,info(i)];
end
interleave=strtrim(interleave);%删除字符串中前后的空格
%读取图像文件
fid = fopen(imagefile, 'r');
data = multibandread(imagefile ,[lines, samples, bands],precision,str2num(headeroffset),interleave,'ieee-le');
% data=double(data);%将影像uint16格式转化为double格式便于显示
% R=data(:,:,1);
% R=data(:,:,2);
% R=data(:,:,3);
% figure
% imagesc(R);%uint8对应double[0,1] uint16对应于double并不是[0,1]的范围，所以显示出来就是一片白
% data=double(data);
% imshow(uint8(data(:,:,1)));
fclose(fid);
end
去极化线性拉伸函数
在这里插入代码片function MAT=LineExtension(MAT)
[row,col] = size(MAT);
%保存每个DN值出现的频率
grayscale = 2^Getgrayscale(MAT)+100;
cum=zeros(1,grayscale);
minDN=0;
%PCA变换后的值有负值
for i=1:row
for j=1:col
if minDN>MAT(i,j)
minDN=MAT(i,j);
end
end
end
%此时cum(1,minDN+abs(minDN)+1)就对应的最小值----实际的cum(1,1)
for i=1:row
for j=1:col
index = round(MAT(i,j));
index=round(index+abs(minDN)+1);
% if index<=grayscale
try
cum(1,index)=cum(1,index)+1;
catch
disp(index);
disp(MAT(i,j));
end
% else
% cum(1,grayscale)=cum(1,grayscale)+1;
% end
end
end
cum=double(cum)/double(row)/col;
final = zeros(1,grayscale);
final(1,1)=cum(1,1);
min=0;
minP=10000;
max=0;
maxP=10000;
for i=2:grayscale
final(1,i)=final(1,i-1)+cum(1,i);
temp=abs(final(1,i)-0.01);
if minP>temp
min=i-1;
minP=temp;
end
temp1=abs(final(1,i)-0.99);
if maxP>temp1
max=i-1;
maxP=temp1;
end
end
%求解得到累积概率之后就进行%2到%98的拉伸------max-abs(minDN),min-abs(minDN)就是2%-98%
for i=1:row
for j=1:col
MAT(i,j)=round(255/(max-min)*(MAT(i,j)-(min-minDN)));
if MAT(i,j)>255
MAT(i,j)=255;
end
if MAT(i,j)<0
MAT(i,j)=0;
end
end
end
% imshow(uint8(MAT));
end
%获取影像的灰度级数
function i = Getgrayscale(A)
max=0;
[row,col]=size(A);
for i=1:row
for j=1:col
if max<A(i,j)
max=A(i,j);
end
end
end
% disp(strcat('最大DN值为',int2str(max)));
for i=8:16
if max<2^i
break
end
end
% disp(strcat('对应的灰度级为',int2str(i)));
end
% function [MIN,MAX]=GetMAXMIN(MAT)
% [row,col]=size(MAT);
% MAX=0
% MIN=10000
% for i=1:row
% for j=1:col
% if MAT(i,j)>MAX
% MAX=MAT(i,j)
% end
% if MAT(i,j)<MIN
% MIN=MAT(i,j);
% end
% end
% end
% end
调用函数的主函数
function test1()
%%
% %Brovey变换
% msidata = read_ENVIimage('Landsat8_OLI_multi.dat');%获取多光谱影像数据，测试数据为4个波段，bsq的存储方式，分辨率为1024*1024
% pandata=read_ENVIimage2('Landsat8_OLI_pan.dat');%获取全色波段的图像矩阵，大概为4096*4096，相对于分辨率提高了4倍
% Out=Brovey(msidata,pandata);
% figure
% for i=1:3
% msidata(:,:,i)=LineExtension(msidata(:,:,i));
% end
% imshow(uint8(msidata(:,:,1:3)));title('msi')%uint16转换到0-1
% figure
% pandata=LineExtension(pandata);
% imshow(uint8(pandata));title('pan');%uint16转换到0-1
% figure
% for i=1:3
% Out(:,:,i)=LineExtension(Out(:,:,i));
% end
% imshow(uint8(Out(:,:,1:3)));title('Brovey');
%%
% %PCA融合变换
% msidata = read_ENVIimage('Landsat8_OLI_multi.dat');%获取多光谱影像数据，测试数据为4个波段，bsq的存储方式，分辨率为1024*1024
% pandata=read_ENVIimage2('Landsat8_OLI_pan.dat');%获取全色波段的图像矩阵，大概为4096*4096，相对于分辨率提高了4倍
% Out=KLFusion(msidata,pandata);
% imshow(uint8(Out));title('PCA');
%%
% % HSI融合变换
% msidata = read_ENVIimage('Landsat8_OLI_multi.dat');%获取多光谱影像数据，测试数据为4个波段，bsq的存储方式，分辨率为1024*1024
% pandata=read_ENVIimage2('Landsat8_OLI_pan.dat');%获取全色波段的图像矩阵，大概为4096*4096，相对于分辨率提高了4倍
% HSIFusion(msidata,pandata);
%%
%Kmeans算法
data = read_ENVIimage3('Landsat8_OLI_multi_classify.dat');
data=double(data);
[~,~,bands]=size(data);
for i=1:bands
data(:,:,i)=LineExtension(data(:,:,i));
end
kinds=input('Please input the kinds that you''d like to classify:');
Kmeans(data,kinds);
figure
imshow(uint8(data(:,:,1:3)));title('origin image');
end
马氏距离与欧式距离
马氏距离(Mahalanobis Distance)是度量学习中一种常用的距离指标，同欧氏距离、曼哈顿距离、汉明距离等一样被用作评定数据之间的相似度指标。但却可以应对高维线性分布的数据中各维度间非独立同分布的问题。
马氏距离(Mahalanobis Distance)是一种距离的度量，可以看作是欧氏距离的一种修正，修正了欧式距离中各个维度尺度不一致且相关的问题。
实际上，马氏距离非常类似于欧式距离，但是呢，中间过了一个方差、协方差的逆阵问题。根据误差传播原理与误差理论可知，方差协方差的逆阵实际上就是对应观测值的权阵。将各个观测值与集群中心之间的差距在各个维度进行归一化，乘以不同的权值之后就获得了更好的归一化之后的马氏距离。相对来说，马氏距离更能代表广义上的距离。具体关于马氏距离可以参考这篇知乎文章，马氏距离在这篇文章中的kmeans并没有使用到，使用的是欧式距离，当然你可稍微改改就可。
851
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
