数字图像处理程序合集上（matlab实现）

---

前言

数字图像处理是测绘科学与技术的基础专业课之一，本篇文章包含了数字图像处理这门课程中常用的功能实现代码合集上集，这个合集将会分为上中下三个部分，当然还有许多不足之处，请大家批判性借鉴。

1.图像读取(读取ENVI软件存储的“.img”、“.dat”格式数据)

path1=input("请输入头文件的路径");
path2=input("请输入数据文件的路径");
data=read_ENVI(path1,path2);
mean=cal_mean(data);
center=cal_center(data);
hist=cal_hist(data);
maxnum=cal_maxnum(data);
cov=cal_cov(data);
r=cal_r(data);
disp(["第",num2str(num1),"波段的像素平均数为:",num2str(av1),",中位数为:",num2str(ca1),",众数为:",num2str(zh1),",方差为:",num2str(fca)])
disp(["第",num2str(num2),"波段的像素平均数为:",num2str(av2),",中位数为:",num2str(ca2),",众数为:",num2str(zh2),",方差为:",num2str(fcb)])
disp(["两个波段像素协方差为:",num2str(xc),",相关系数为:",num2str(xg)])
function data=read_ENVI(hdrfilename,imgfilename)
fid = fopen(hdrfilename, 'r');
info = fread(fid,'char=>char');
info=info';%默认读入列向量，须要转置为行向量才适于显示
fclose(fid);
%查找列数
a=strfind(info,'samples = ');
b=length('samples = ');
c=strfind(info,'lines');
samples=[];
for i=a+b:c-1
    samples=[samples,info(i)];
end
samples=str2num(samples);
%查找行数
a=strfind(info,'lines   = ');
b=length('lines   = ');
c=strfind(info,'bands');
lines=[];
for i=a+b:c-1
    lines=[lines,info(i)];
end
lines=str2num(lines);
%查找波段数
a=strfind(info,'bands   = ');
b=length('bands   = ');
c=strfind(info,'header offset');
bands=[];
for i=a+b:c-1
    bands=[bands,info(i)];
end
bands=str2num(bands);
%查找数据类型
a=strfind(info,'data type = ');
b=length('data type = ');
c=strfind(info,'interleave');
datatype=[];
for i=a+b:c-1
    datatype=[datatype,info(i)];
end
datatype=str2num(datatype);
precision=[];
switch datatype
    case 1
        precision='uint8=>uint8';%头文件中datatype=1对应ENVI中数据类型为Byte，对应MATLAB中数据类型为uint8
    case 2
        precision='int16=>int16';%头文件中datatype=2对应ENVI中数据类型为Integer，对应MATLAB中数据类型为int16
    case 12
        precision='uint16=>uint16';%头文件中datatype=12对应ENVI中数据类型为Unsighed Int，对应MATLAB中数据类型为uint16
    case 3
        precision='int32=>int32';%头文件中datatype=3对应ENVI中数据类型为Long Integer，对应MATLAB中数据类型为int32
    case 13
        precision='uint32=>uint32';%头文件中datatype=13对应ENVI中数据类型为Unsighed Long，对应MATLAB中数据类型为uint32
    case 4
        precision='float32=>float32';%头文件中datatype=4对应ENVI中数据类型为Floating Point，对应MATLAB中数据类型为float32
    case 5
        precision='double=>double';%头文件中datatype=5对应ENVI中数据类型为Double Precision，对应MATLAB中数据类型为double
    otherwise
        error('invalid datatype');%除以上几种常见数据类型之外的数据类型视为无效的数据类型
end
%查找数据格式
a=strfind(info,'interleave = ');
b=length('interleave = ');
c=strfind(info,'sensor type');
interleave=[];
for i=a+b:c-1
    interleave=[interleave,info(i)];
end
    interleave=strtrim(interleave);%删除字符串中的空格
%读取图像文件
fid = fopen(imgfilename, 'r');
data = multibandread(imgfilename ,[lines, samples, bands],precision,0,interleave,'ieee-le');
data= double(data);
end

2.统计与描述

2.1 均值

function mean=cal_mean(A) % 计算矩阵，给的是三维矩阵
[row,col,band]=size(A);
A=reshape(A,[row*col,band]);
sum=zeros(1,band);
for i=1:band
    for j=1:row*col
        sum(1,i)=sum(1,i)+A(j,i);
    end
end
mean=sum./(col*row);
end

2.2中值加上排序

function [A,center]=cal_center(A,sort) % sort为排序的方式，1代表正序，2代表倒序
[row,col,band]=size(A);
A=reshape(A,[row*col,band]);
samples=row*col;
if sort==1
    for i=1:band
     for j=1:samples
        for k=1:samples-j
            if A(k,i)>A(k+1,i)
                temp=A(k,i);
                A(k,i)=A(k+1,i);
                A(k+1,i)=temp;
            end
        end
     end
    end % sort=1是从小到大排序
elif sort==2
    for i=1:bnad
     for j=1:samples
        for k=1:samples-j
            if A(k,i)<A(k+1,i)
                temp=A(k,i);
                A(k,i)=A(k+1,i);
                A(k+1,i)=temp;
            end
        end
     end
   end
end % sort=2是从大到小排序

for i=1:band
    if mod(samples,2)==0
        center=(A(samples/2,i)+A(samples/2+1,i))/2;
    else
        center=A(samples/2+1,i);
    end % 根据samples是奇数还是偶数，找到最终的那个数作为中数
end
end

2.3 众数

function [index,maxnum]=cal_maxnum(A) % 求解众数，同时输出出现的次数
hists=cal_hist(A);% 调用了自己编写的cal_hist函数得到频数分布表
hists=hists(:,2,:);
[row,~,band]=size(hists);
maxnum=zeros(1,band);
for i=1:band
    for j=1:row
        if hists(j+1,1,i)>hists(j,1,i)
            maxnum=j+1;
            index=hists(j+1,1,i);
        end
    end
end
end

2.4 方差与协方差矩阵

function cov=cal_cov(A) %输出的是方差-协方差矩阵，主对角线上元素是方差，cov(i,j)表示第i个波段与第j个波段的协方差
[row,col,band]=size(A);
samples=row*col;
mean=repmat(cal_mean(A),[samples,1]);
new=zeros(samples,band);
for i=1:band
    B=reshape(A(:,:,i),[samples,1]);
    new(:,i)=B;
end
cov=zeros(band,band);
for i=1:band
    for j=1:band
        sum1=new(:,i)-mean(:,i);
        sum2=new(:,j)-mean(:,j);
        sum3=sum(sum1.*sum2);
        cov(i,j)=sum3./samples;
    end
end
end

2.5 相关系数

function r=cal_r(A) % 根据协方差-方差矩阵得到各个波段数据的相关系数，主对角线的值为1
cov=cal_cov(A);
[row,col]=size(cov);
r=zeros(row,col);
for i=1: row
    for j=1:col
        r(i,j)=cov(i,j)/(pow2(cov(i,i),0.5)*pow2(cov(j,j),0.5));
    end
end
end

2.6 直方图/累计直方图

function hists=cal_hist(A) %返回的是一个256*2*波段数的矩阵，第一列是分布概率，第二列是累计分布概率
[row,col,band]=size(A);
B=A;
samples=row*col;
A=reshape(A,[samples,band]);
maxs=max(A);
hists=zeros(256,2,band);
for i=1:band
    for j=1:samples
        for k=1:maxs(i)
            if A(j,i)==k-1
                hists(k,1,i)=hists(k,1,i)+1;
            end
        end
    end
end
hists=hists./samples;
hists(1,2,:)=hists(1,1,:);
for i=1:255
    hists(i+1,2,:)=hists(i,2,:)+hists(i+1,1,:);
end
hudu=0:255;
for i=1:band % 分波段画出直方图
    figure(i)
    histogram(B(:,:,i))
end
for i=1:band % 分波段画出累计直方图
    figure(i+3)
    bar(hudu,hists(:,2,i))
end
end