Matlab 高光谱影像信息熵/信噪比计算

基于matlab实现了高光谱图像波段信噪比和信息熵的计算

1. 文件导入：图像格式使用ENVI导出的img hdr格式，参考Matlab实现高光谱读取进行的修改。时间伧俗，程序适用性不强，影像的参数需要用记事本打开hdr手动修改程序的文件。 —— "Float"修改成你的img例如，文件名example.img —— bands 、lines 、samples 、interleave可以从hdr文件里找到 —— ‘float32=>float32’ 对应hdr 里的data type =4.代表不同数据格式的数据，如果这里的参数不同，需要修改，可以参考这里switch datatype case位置的代码 对应修改。

2. 波段信息熵计算：当导入数据时float因此，将其分为N组离散统计组概率，然后应用公式计算波段信息熵。 注：imhist函数对于doubel类型分组的范围在0-1之间，因此数据值需要在此范围内。

3. 波段信噪比计算：主要参考：ENVI IDL实现波段信噪比计算matlab版本，详细原理也可以参考本文

-------------有问题可以互相讨论和交流--------------

%% 文件导入 tic fprintf("正在打开。。。"); 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:bands     tic     temp = 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 % 遥感图像计算信噪比的方法不同 % 参考资料: % https://blog.csdn.net/qq_33356563/article/details/82081081 % https://www.docin.com/p-1473544620.html fprintf("计算信噪比。。。。"); width=4; SNR = []; for k = 1:bands     tic     temp = data(:,:,k);     % 1.边缘提取-基于Canny结果是二值图像     mask=edge(temp,'canny');%[0.4,0.6],0.6     % 2.去除边缘块-按规定字块尺寸(4)x4)将整个图像分块，     %   统计每个子块是否包含边缘值，如果有，将子块拆除，不再参与以下信噪比估算。     %逐块计算标准差和均值     nosie_subset =zeros(uint16(lines/width),uint16(samples/width));     singal_subset = zeros(uint16(lines/width),uint16(samples/width));     for i=0:uint16(lines/width)-2         for 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);             end         end     end     % 3.局部方差法估计噪声值     singal_subset = singal_subset(singal_subset~=0);     nosie_subset = nosie_subset(nosie_subset~=0);     % 在局部标准差最小值与平均值的1.150个区间分为两倍，     % 各子块按标准差大小落入相应区间，计算直方图。     binstart = min(nosie_subset);     binend = mean(nosie_subset)*1.2;      binsize =(binend-binstart)/150;     bincount = [];     for i=1:150         bin = 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~=)；%信号估计值     sn = 10*log(singal/noise);% 单位dB          fprintf('第%d个波段耗时%f s\n',k,toc)     SNR = [SNR ,sn]; end %可视化 bar(SNR) 

https://blog.csdn.net/qq_33356563/article/details/82081081 https://www.docin.com/p-1473544620.html