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.youkuaiyun.com/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.边缘块剔除--------按照规定字块尺寸（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)-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.2倍之间划分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~=0));%信号估计值
    sn = 10*log(singal/noise);% 单位dB
    
    fprintf('第%d个波段用时%f s\n',k,toc)
    SNR = [SNR ,sn];
end
%可视化
bar(SNR)

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