clear;close all;clc;
% 读取高光谱图像
X = multibandread('F:/baker1-5_ref.raw',...
     [1133 640 512],'float32',0,'bil','ieee-le'); %[lines samples bands]
[r,c,n] = size(X);
% X=mat2gray(X,[0,1800]);
X_oneband = X(:,:,127);% 显示512幅图中的第127band光谱图
str = '第127band光谱图';
%figure,imshow(X_oneband),title(str);
% 确定阈值范围
figure,imhist(X_oneband)
% 去除阈值大于0.25小于等于0.15的背景
X_oneband(X_oneband > 0.90) = 0;
X_oneband(X_oneband <= 0.3)= 0;
X_oneband(X_oneband > 0.3 & X_oneband < 0.90) = 1;
% Bw = im2bw(X_oneband);
Bw = X_oneband;
%figure,imshow(Bw),title('二值图像');
% 联通区域标记 
[L,num] = bwlabel(Bw,8);
stats = regionprops(L);
% 求取所有联通区域
Ar = cat(1,stats.Area);
thres(r)=30;
thres(c)=30;
ind = find(Ar > thres);
Bw_circle = ismember(L,ind);
%figure,imshow(Bw_circle),title('所有联通区域');
% 填充图像
Bw_circleall = imfill(Bw_circle,'holes');
figure,imshow(Bw_circleall),title('填充图像');
% 形态学操作
se = strel('disk',6);
Bw_erode = imerode(Bw_circleall,se);
%figure,imshow(Bw_erode),title('腐蚀图像');

% 读取平均数据
Re1 = reshape(Bw_erode,r*c,1);% 掩膜板
Re2 = reshape(X,r*c,n);% 原图像
%Re2 = mat2gray(Re2);% 归一化
PSF=fspecial('average',12);      %创建一个结构体，滤波的模板大小为12*12
J=imfilter(Re2,PSF);%均值滤波
% 平均数据
[p,m]=find(Re1==1); 
spect=mean(J(p,:));
% 画光谱图
specter = spect(:,8:486);
bandband = 899.83:8019/4780:1701.73;
%figure(7),plot(bandband,specter),title('掺假光谱');

disp(specter);
a=specter(:);
disp(a');
xlswrite('G:\matlab2023.xlsx',a)