Matlab遥感数字图像处理——统计值计算

Matlab遥感数字图像处理——统计值计算

clear,clc;   %清空工作区，清屏
result = IMG_Read('NVIS.dat');
% result = IMG_Read('Sandiego.img');
% result = IMG_Read('TM.dat');
% result = IMG_Read('can.dat');
% result = IMG_Read('TM-30m.dat');
% result = IMG_Read('Beijing.dat');

datatype = cell2mat(result(1));
Image = cell2mat(result(2));
[lines,samples,bands] = size(Image);

IMG = Image(:,:,1);        %只读取图像的一个波段
IMAGE = reshape(IMG,1,[]); %将二维矩阵转换为一维数组
IMGnum = samples*lines;

IMGsum = 0;
IMGVAR = 0;
IMGCOV = 0;

fid=fopen('统计结果.txt','w');

%求单波段数据的均值
IMGmean = mean(IMAGE)
% for i=1:IMGnum
%     IMGsum = IMGsum + IMAGE(i);
% end
% IMGmean = IMGsum/IMGnum;
fprintf(fid,'%7s %12s\r\n','统计量','统计值');
fprintf(fid,'%7s %17f\r\n','均值',IMGmean);

%对单波段数据排序
IMAGE1 = IMAGE;
IMAGE1 = sort(IMAGE1);
% for j = 1:IMGnum
%  for i = 1:IMGnum-j
%         if(IMAGE1(i)>IMAGE1(i+1))
%             temp = IMAGE1(i);
%             IMAGE1(i) = IMAGE1(i+1);
%             IMAGE1(i+1) = temp;
%         end
%     end
% end

% 求单波段数据的中值
IMGmedian = median(IMAGE1)
% if mod(IMGnum,2)==0
%     IMGmedian = (IMAGE1(IMGnum/2)+IMAGE1(IMGnum/2+1))/2;
% else 
%     IMGmedian = IMAGE1(fix(IMGnum/2));
% end
fprintf(fid,'%7s %17f\r\n','中值',IMGmedian);

%直方图统计
IMGmax = max(IMAGE1)
numbers = zeros(2,IMGmax);
for k = 1:IMGmax+1
    numbers(1,k) = k-1;  
    numbers(2,k)= 0;
end
for i = 1:samples*lines
        for k = 1:IMGmax+1
             if IMAGE1(1,i)== numbers(1,k)
                 numbers(2,k) = numbers(2,k)+1;
             end
        end
end

%累计直方图统计
NUMBER = zeros(2,fix(IMGmax));
NUMBER(2,1) = numbers(2,1);
for k = 2:IMGmax+1
    NUMBER(1,k) = k-1;  
    NUMBER(2,k) = NUMBER(2,k-1) + numbers(2,k); 
end

% 求单波段数据的众数
IMGmode = mode(IMAGE1(1,:))
fprintf(fid,'%7s %17f\r\n','众数',IMGmode);
% for i = 1:IMGmax
%     if(numbers(2,i)== max(numbers(2,:)))
%         IMGmode = i-1;
%     end
% end

% 求单波段数据的方差
IMGvar = var(IMAGE)
% for i = 1:IMGnum
%    IMGVAR = IMGVAR+ (IMAGE(i)-IMGmean) ^2;
% end
% IMGvar = IMGVAR/IMGnum
fprintf(fid,'%7s %17f\r\n','方差',IMGvar);

%前两个波段的协方差
IMG2 = Image(:,:,2);
IMAGE2 = reshape(IMG2,1,[]);
IMGcov = cov(IMAGE,IMAGE2)
% IMGmean2 = mean(IMAGE2);
% for i = 1:IMGnum
%    IMGCOV = IMGCOV+ (IMAGE(i)-IMGmean) * (IMAGE2(i)-IMGmean2);
% end
% IMGcov = IMGCOV/IMGnum
fprintf(fid,'%7s\r %16f %10f\r\n %26f %10f\r\n','协方差',IMGcov);

%前两个波段的相关系数
IMGcorrcoef = corrcoef(IMG,IMG2)
% TMGcorrcoef = IMGcov/((IMGvar^0.5)*(var(IMAGE2))^0.5)
fprintf(fid,'%7s\r %14f %10f\r\n %25f %10f\r\n','相关系数',IMGcorrcoef);
fclose(fid);

%绘制单波段数据的直方图
figure(2)
subplot(1,3,1);
hist(IMAGE1,IMGmax)
title('第一个波段灰度直方图');
% IMGhist = hist(IMAGE,IMGmax);

subplot(1,3,2);
hist(IMAGE2,max(IMAGE2));
title('第二个波段灰度直方图');
% IMGhist2 = hist(IMAGE2,max(IMAGE2));

subplot(1,3,3);
IMG3 = Image(:,:,3);
IMAGE3 = reshape(IMG3,1,[]);
hist(IMAGE3,max(IMAGE2))
title('第三个波段灰度直方图');
% IMGhist3 = hist(IMAGE3,max(IMAGE2));