高光谱影像信息熵/信噪比计算
基于matlab实现了高光谱影像个波段信噪比和信息熵的计算
文件导入：影像格式使用的ENVI导出的img+hdr格式，参考的Matlab实现高光谱读取进行的修改。时间伧俗，程序适用性不强，影像的参数需要用记事本打开hdr文件，手动修改程序。
—— "Float"修改成你的img文件名，例如example.img
—— bands 、lines 、samples 、interleave可以从hdr文件里找到
—— ‘float32=>float32’ 对应hdr 里的data type =4。代表数据不同的数据格式，此处参数不同的话需修改，可以参考此处switch datatype case位置的代码 进行对应修改。
波段信息熵计算：因为导入的数据时float型，所以将其分为N组离散化统计组概率，进而才可以套用公式计算波段信息熵。
注：imhist函数对于doubel类型分组的范围为0-1之间，所以需要数据的值在这个范围内。
波段信噪比计算：主要参考ENVI IDL波段信噪比计算实现的matlab版本，详细的原理还可参考这篇论文
-------------有问题可互相探讨交流--------------
程序
%% 文件导入ticfprintf("正在打开。。。");fid = fopen("Float", 'r');bands = 274;lines = 2437;samples = 2379;interleave = 'bsq';data = multibandread("Float" ,[lines, samples, bands],'float32=>float32',0,interleave,'ieee-le');data= double(data);fprintf('文件导入用时%f s\n', toc)%% 显示某波段imshow(data(:,:,100));%% 计算信息熵% dataMin = min(min(min(data)));% dataMax = max(max(max(data)));N = 2^14; %直方分组数量Entropy = [];for i = 1:bandstictemp = data(:,:,i);% 统计单波段直方图，并计算每个区间内的出现概率，0区间是背景不参与统计[counts,x] = imhist(temp , N);P = counts(2:end)/sum(counts(2:end));% 使用公式进行计算，L2P = log2(P);L2P(L2P==-inf)=0;H = -sum(P.*L2P);fprintf('第%d个波段用时%f s\n',i,toc)Entropy = [Entropy,H];end%归一化到0-1之间Entropy = Entropy/max(Entropy);%可视化bar(Entropy)%% 计算信噪比SNR ---基于局部方差法计算遥感影像的信噪比 使用ENVI IDL% 遥感影像计算信噪比的方法不太一样% 参考资料:% /qq_33356563/article/details/82081081% /p-1473544620.htmlfprintf("正在计算信噪比。。。");width=4;SNR = [];for k = 1:bandstictemp = data(:,:,k);% 1.边缘提取--------基于Canny算子对图像进行边缘提取，结果为二值图像mask=edge(temp,'canny');%[0.4,0.6],0.6% 2.边缘块剔除--------按照规定字块尺寸（4x4）对整个图像进行分块，% 统计每一个子块中是否包含有边缘值，若有则将该子块剔除，不再参加后面的信噪比估算。%逐块计算标准差和均值nosie_subset =zeros(uint16(lines/width),uint16(samples/width));singal_subset = zeros(uint16(lines/width),uint16(samples/width));for i=0:uint16(lines/width)-2for j=0:uint16(samples/width)-2%如果当前块含有边缘，则进入下一像元tmask=mask(i*width+1:(i+1)*width,j*width+1:(j+1)*width);if sum(sum(tmask)) ==0 %当前分块数据的标准差及均值tdata=temp(i*width+1:(i+1)*width,j*width+1:(j+1)*width);nosie_subset(i+1,j+1) = std2(tdata);singal_subset(i+1,j+1) = mean2(tdata);endendend% 3.局部方差法估算噪声值singal_subset = singal_subset(singal_subset~=0);nosie_subset = nosie_subset(nosie_subset~=0);% 在局部标准差最小值与平均值的1.2倍之间划分150个区间，% 按标准差大小将各子块落入到相应的区间，以此计算得到直方图。binstart = min(nosie_subset);binend = mean(nosie_subset)*1.2; binsize =(binend-binstart)/150;bincount = [];for i=1:150bin = sum((nosie_subset>=binstart+binsize*(i-1))&(nosie_subset<binstart+binsize*i));bincount = [bincount,bin]; end% 根据直方图统计出包含子块最多的区间，计算区间内标准差的平均值作为噪声估计值。 [~,i]= max(bincount);noise = mean(nosie_subset((nosie_subset>=binstart+binsize*(i-1))&(nosie_subset<binstart+binsize*i)));% 4.信噪比计算--------统计剔除边缘块后的像元平均值作为估计值；波段平均标准差除以该波段的平均值，得到信噪比singal = mean(singal_subset(singal_subset~=0));%信号估计值sn = 10*log(singal/noise);% 单位dBfprintf('第%d个波段用时%f s\n',k,toc)SNR = [SNR ,sn];end%可视化bar(SNR)
/qq_33356563/article/details/82081081
/p-1473544620.html
