| Matlab遥感数字图像处理函数汇总 | 您所在的位置：网站首页 › 利用遥感影像如何进行制图工作 › Matlab遥感数字图像处理函数汇总 |
|
Matlab遥感数字图像处理函数汇总
遥感数字图像处理 是通过计算机图像处理系统 对遥感图像中的像素进行的系列操作过程。 （以下函数均针对每一个波段分开进行处理，如需综合处理，用reshape函数将图像矩阵转换成 向量(行×列×波段数) 即可） 温馨提示：‘文章有点长，可按目录定位！’ 文章目录 Matlab遥感数字图像处理函数汇总1 图像读取1.1 读取遥感影像hdr格式头文件1.2 读取img格式或dat格式的遥感图像 2 图像统计与描述2.1 图像均值函数2.2 图像中值函数2.3 图像累计直方图函数2.4 图像众数函数2.5 图像的协方差矩阵2.6 图像的相关系数矩阵 3 图像增强处理3.1 直方图线性拉伸3.2 直方图均衡化3.3 直方图匹配3.4 中值滤波3.5 均值滤波3.6 Sobel算子锐化3.7 Prewitt算子锐化3.8 RGB to HSI3.9 HSI to RGB3.10 PCA变换 4 影像融合处理4.1 基于Brovey变换的影像融合4.2 基于乘积变换的影像融合4.3基于PCA的影像融合4.4 基于HSI变换的影像融合 5 影像特征提取5.1 迭代阈值影像分割算法5.2 区域生长影像分割算法 6 影像分类6.1 K-Means聚类分割6.2 最短距离法分类6.3 最大似然法分类 1 图像读取 1.1 读取遥感影像hdr格式头文件 function parameter =HDR_Read(HDRfilename) %设置函数HDR_Read，返回值为parameter，参数为HDRfilename fileID = fopen(HDRfilename, 'r'); %以只读方式打开.hdr头文件 %一直读取头文件的内容直至文件末尾 while ~feof(fileID) fileline=fgetl(fileID); C = cellstr(strsplit(fileline,'= ')); %按"="分割每一行 if(C(1)=="samples ") samples=str2double(cell2mat(C(2))); %存储列数samples end if(C(1)=="lines ") lines=str2double(cell2mat(C(2))); %存储行数lines end if(C(1)=="bands ") %存储波段数bands bands=str2double(cell2mat(C(2))); end if(C(1)=="interleave ") %存储envi中数据存储方式 interleave=cell2mat(C(2)); end if(C(1)=="data type ") %存储数据类型 datatype=cell2mat(C(2)); end end fclose(fileID); %关闭文件的读取 %定义ENVI中的数据类型格式 data_type = {'uint8' 'int16' 'int32' 'float32' 'float64' 'uint16' 'uint32' 'uint64'}; %将envi中不同的数据格式转换成matlab中的特定格式 switch datatype case '1' datatype = cell2mat(data_type(1)); %灰度范围0-255 case '2' datatype = cell2mat(data_type(2)); %16位整数 case '3' datatype = cell2mat(data_type(3)); %32位整数 case '4' datatype = cell2mat(data_type(4)); %32位浮点数 case '5' datatype = cell2mat(data_type(5)); %64位浮点数 case '6' datatype = cell2mat(data_type(6)); %灰度范围0-2^16 case '7' datatype = cell2mat(data_type(7)); %灰度范围0-2^32 case '8' datatype = cell2mat(data_type(8)); %灰度范围0-2^64 end %将各数据存入返回值参数parameter中 parameter = {samples,lines,bands,interleave,datatype}; end 1.2 读取img格式或dat格式的遥感图像 function Re=IMG_Read(IMGfilename) HDRfilename =strcat(strtok(IMGfilename,'.'),'.hdr'); %得到相同名字的头文件 parameter = HDR_Read(HDRfilename); %读取头文件中的参数信息 samples = cell2mat(parameter(1)); %得到列数 lines = cell2mat(parameter(2)); %得到行数 bands = cell2mat(parameter(3)); %获得波段数 interleave = cell2mat(parameter(4)); %获得图像存储格式 datatype = cell2mat(parameter(5)); %获取数据格式 fid = fopen(IMGfilename, 'r'); %以只读方式打开img文件 Image=fread(fid,datatype); %按指定格式读取img图像文件 %针对不同的图像存储格式，处理的到图像矩阵 switch interleave case 'bsq' %当存储格式为bsq时,按波段存储 Image=reshape(Image,[samples,lines,bands]);%将图像读入三维矩阵 for i=1:bands IMG(:,:,i)=Image(:,:,i)'; end Image=IMG; case 'bil' %当存储格式为bil时，按行存储 IMG=zeros(lines,samples,bands); count=1; for row =1:lines for i=1:bands IMG(row,:,i) =Image((count-1)*samples+1:count*samples); count=count+1; end end Image=IMG; case 'bip' %当存储格式为bip时，按像元存储 IMG=zeros(lines,samples,bands); count=1; for row=1:lines for col=1:samples IMG(row,col,:)=Image((count-1)*bands+1:count*bands); count=count+1; end end Image=IMG; otherwise error('格式错误，无法读取图像！'); end Re = {samples,lines,bands,datatype,Image}; fclose(fid); subplot(2,2,1); imgmul=cat(3,Image(:,:,3),Image(:,:,2),Image(:,:,1));%合成3维矩阵 colormap(colorcube); imshow(uint8(imgmul)); colorbar; title('前三个波段组合图像'); subplot(2,2,2); img1=Image(:,:,1); colormap(gray); imshow(uint8(img1)); colorbar; title('第一个波段灰度图'); subplot(2,2,3); img2=Image(:,:,2); colormap(gray); imshow(uint8(img2)) colorbar; title('第二个波段灰度图'); subplot(2,2,4); img3=Image(:,:,3); colormap(gray); imshow(uint8(img3)) colorbar; title('第三个波段灰度图'); % axis off; 2 图像统计与描述 2.1 图像均值函数 %函数名称为Image_Mean,输入参数Image,输出参数Mean function [Mean] = Image_Mean(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %将三维矩阵转换成二维矩阵，方便计算 Image1 = reshape(Image,[m*n,bands]); %初始化求和矩阵 Sum = zeros(bands,1); %计算每个波段所有像元的灰度之和 for k = 1:bands for i = 1:m*n Sum(k) = Sum(k) + Image1(i,k); end end %计算每个波段的灰度平均值 Mean = Sum /(m*n); end 2.2 图像中值函数 %函数名称为Image_Median,输入参数Image,输出参数Median function [Median] = Image_Median(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %将三维矩阵转换成二维矩阵，方便计算 Image1 = reshape(Image,[m*n,bands]); %对每一个波段进行冒泡排序 for k = 1:bands for i = 1:m*n for j = 1:m*n-i if(Image1(j,k)>Image1(j+1,k)) temp = Image1(j,k); Image1(j,k) = Image1(j+1,k); Image1(j+1,k) = temp; end end end end %去排序后每个波段中间值作为中值 for k = 1:bands if(mod(m*n,2) == 0) Median(k) = Image1(m*n/2,k)/2 + Image1(m*n/2+1,k)/2; else Median(k) = Image1(m*n/2+1,k); end end end 2.3 图像累计直方图函数 %函数名称为Image_Hist,输入参数Image,输出参数Hist function [Hist] = Image_Hist(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %将三维矩阵转换成二维矩阵，方便计算 Image1 = reshape(Image,[m*n,bands]); %初始化三维矩阵，行表示256种灰度，列表示灰度值、个数、累计个数 Hist = zeros(256,2,bands); %求每个波段中每个灰度值的个数 for k = 1:bands for i = 1:256 for j = 1:m*n if(Image1(j,k) == i-1) Hist(i,1,k) = Hist(i,1,k) + 1; end end end end %转换为频率直方图 Hist = Hist./(m*n); %求每个波段每个灰度的累计个数 Hist(1,2,:) = Hist(1,1,:); for k = 1:bands for i = 2:256 Hist(i,2,k) = Hist(i-1,2,k) + Hist(i,1,k); end end end 2.4 图像众数函数 %函数名称为Image_Mode,输入参数Image,输出参数Mode function [Mode] = Image_Mode(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %将三维矩阵转换成二维矩阵，方便计算 Image1 = reshape(Image,[m*n,bands]); %计算直方图矩阵 Hist = Image_Hist(Image); %另初始最大个数为0 Max = 0 ; %初始众数矩阵，每一行是一个波段的众数 Mode = zeros(bands,1); %计算每个波段灰度值的众数 for k = 1:bands for i = 1:256 if(Hist(i,1,k)>Max) Max = Hist(i,1,k); Mode(k) = i-1; end end end end 2.5 图像的协方差矩阵 %函数名称为Image_Cov,输入参数Image,输出参数Cov function [Cov] = Image_Cov(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %将三维矩阵转换成二维矩阵，方便计算 %只有三个波段的情况下 Image1 = reshape(Image(:,:,1),[m*n,1]); Image2 = reshape(Image(:,:,2),[m*n,1]); Image3 = reshape(Image(:,:,3),[m*n,1]); IM = [Image1;Image2;Image3]; %求每一个波段灰度值的均值 Mean = Image_Mean(Image); %初始化协方差矩阵 Cov = zeros(bands); %计算协方差矩阵 for i = 1:bands for j = 1:bands for k = 1:m*n Cov(i,j) = Cov(i,j)+(IM(i,k)-Mean(i))*(IM(j,k)-Mean(k)) ; end Cov(i,j) = Cov(i,j); end end end 2.6 图像的相关系数矩阵 %函数名称为Image_R,输入参数Image,输出参数R function [R] = Image_R(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %计算协方差矩阵 Cov = Image_Cov(Image); %初始化相关系数矩阵 R = zeros(bands); %计算相关系数 for i = 1:bands for j = 1:bands R = Cov(i,j)/(sqrt(Cov(i,i))*sqrt(Cov(j,j))); end end end 3 图像增强处理 3.1 直方图线性拉伸 %函数名称为Image_LinearStretch,输入参数Image,输出参数IMAGE function [IMAGE] = Image_LinearStretch(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %初始化计算矩阵 Image1 = reshape(Image,[m*n,bands]); %计算每个波段灰度的最大最小值 for k = 1:bands Max(k,1) = max(max(Image(:,:,k))); Min(k,1) = min(min(Image(:,:,k))); end %线性拉伸 for k = 1:bands for i = 1:m*n Image1(i,k) = (Image1(i,k)-Min(k))*255/(Max(k)-Min(k)) + 0; end end %将二维矩阵转换为三维矩阵 IMAGE = reshape(Image1,[m,n,bands]); %画图，左右分别表示原图和处理后的图像 figure(1); subplot(1,2,1); imshow(uint8(Image)); title('原始图像'); subplot(1,2,2); imshow(uint8(IMAGE)); title('线性拉伸后的图像'); end 3.2 直方图均衡化 %函数名称为Image_HistogramEqualization,输入参数Image,输出参数IMAGE function [IMAGE] = Image_HistogramEqualization(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %初始化计算矩阵 IMAGE = Image; %计算图像的直方图矩阵 Hist = Image_Hist(Image); %建立灰度映射函数 for k = 1:bands for i = 1:256 P(i,k) = round(255* Hist(i,2,k)); end end %计算直方图均衡化后的矩阵 for k = 1:bands for i = 1:m for j = 1:n IMAGE(i,j,k) = P(IMAGE(i,j,k)+1,k); end end end %画图，左右分别表示原图和处理后的图像 figure(1); subplot(1,2,1); imshow(uint8(Image)); title('原始图像'); subplot(1,2,2); imshow(uint8(IMAGE)); title('直方图均衡化后的图像') end 3.3 直方图匹配 %函数名称为Image_HistogramMatching,输入参数Image1，Image2,输出参数IMAGE function [IMAGE] = Image_HistogramMatching(Image1,Image2) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image1); %初始化计算矩阵 IMAGE = Image1; %计算图像的直方图矩阵 Hist1 = Image_Hist(Image1); Hist2 = Image_Hist(Image2); %建立灰度映射函数 for k = 1:bands for i = 1:256 [vv,index] = min(abs(Hist1(i,2,k)-Hist2(:,2,k))); P(i,k) = index -1; end end %计算直方图均衡化后的矩阵 for k = 1:bands for i = 1:m for j = 1:n IMAGE(i,j,k) = P(IMAGE(i,j,k)+1,k); end end end %画图 figure(1); subplot(1,3,1); imshow(uint8(Image1)); title('原始图像'); subplot(1,3,2); imshow(uint8(Image2)); title('参照图像'); subplot(1,3,3); imshow(uint8(IMAGE)); title('直方图匹配化后的图像') end 3.4 中值滤波 %函数名称为Image_MedianFilter,输入参数Image,输出参数IMAGE function [IMAGE] = Image_MedianFilter(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %加高斯噪声 Image1 = imnoise(uint8(Image),'gaussian'); %初始化矩阵 IMAGE = Image1; %定义模板大小，假设模板大小3×3 A = 1; %中值滤波 for k = 1:bands for i = 1+A:m-A for j = 1+A:n-A temp = Image1(i-A:i+A,j-A:j+A,k); IMAGE(i,j,k) = round(median(median(temp))); end end end %画图 figure(1); subplot(1,3,1); imshow(uint8(Image)); title('原始图像'); subplot(1,3,2); imshow(uint8(Image1)); title('加噪声以后的图像'); subplot(1,3,3); imshow(uint8(IMAGE)); title('中值滤波后的图像') end 3.5 均值滤波 %函数名称为Image_Mean,输入参数Image,输出参数IMAGE function [IMAGE] = Image_Mean(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %加高斯噪声 Image1 = imnoise(uint8(Image),'gaussian'); %定义模板大小，假设模板大小3×3 A = 1; %初始化矩阵 IMAGE = Image1; %均值滤波 for k = 1:bands for i = 1+A:m-A for j = 1+A:n-A temp = Image1(i-A:i+A,j-A:j+A,k); IMAGE(i,j,k) = round(mean(mean(temp))); end end end %画图 figure(1); subplot(1,3,1); imshow(uint8(Image)); title('原始图像'); subplot(1,3,2); imshow(uint8(Image1)); title('加噪声以后的图像'); subplot(1,3,3); imshow(uint8(IMAGE)); title('均值滤波后的图像') end 3.6 Sobel算子锐化 %函数名称为Image_Sobel,输入参数Image,输出参数IMAGE function [IMAGE] = Image_Sobel(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %定义模板大小，假设模板大小3×3 A = 1; %定义Sobel算子x，y方向矩阵 Sobelx = [-1 -2 -1;0 0 0;1 2 1]; Sobely = [-1 0 1;-2 0 2;-1 0 1]; %初始化矩阵 Image1 = zeros(m,n,bands); IMAGE = Image; %Sobel算子 for k = 1:bands for i = 1+A:m-A for j = 1+A:n-A temp = Image(i-A:i+A,j-A:j+A,k); Image1(i,j,k) = abs(sum(sum(Sobelx.*temp)))+abs(sum(sum(Sobely.*temp))); end end end IMAGE = Image + Image1; %画图，左右分别表示原图和两幅处理后的图像 figure(1); subplot(1,3,1); imshow(uint8(Image)); title('原始图像'); subplot(1,3,2); imshow(uint8(Image1)); title('边缘提取图像'); subplot(1,3,3); imshow(uint8(IMAGE)); title('Sobel算子锐化后的图像') end 3.7 Prewitt算子锐化 %函数名称为Image_Prewitt,输入参数Image,输出参数IMAGE function [IMAGE] = Image_Prewitt(Image) %获取矩阵的行、列、波段数 [m,n,bands] = size(Image); %定义模板大小，假设模板大小3×3 A = 1; %定义Prewitt算子x，y方向矩阵 Prewittx = [-1 -2 -1;0 0 0;1 2 1]; Prewitty = [-1 0 1;-2 0 2;-1 0 1]; %初始化矩阵 Image1 = zeros(m,n,bands); IMAGE = Image; %Sobel算子 for k = 1:bands for i = 1+A:m-A for j = 1+A:n-A temp = Image(i-A:i+A,j-A:j+A,k); Image1(i,j,k) = abs(sum(sum(Prewittx.*temp)))+abs(sum(sum(Prewitty.*temp))); end end end IMAGE = Image + Image1; %画图，左中右分别表示原图和两幅处理后的图像 figure(1); subplot(1,3,1); imshow(uint8(Image)); title('原始图像'); subplot(1,3,2); imshow(uint8(Image1)); title('边缘提取图像'); subplot(1,3,3); imshow(uint8(IMAGE)); title('Prewitt锐化后的图像') end 3.8 RGB to HSI %函数名称为Image_RGB2HSI,输入参数Image,输出参数HSI function HSI = Image_RGB2HSI(Image) %RGB->HSI变换 [lines,samples,bands] = size(Image); HSI = zeros(lines,samples,3); %用三维向量分别存储HSI（色度，饱和度，亮度） for i = 1:lines for j = 1:samples if( Image(i,j,1)== Image(i,j,2) && Image(i,j,2) ==Image(i,j,3)) JD = 0; else JD = 180*acos((2*Image(i,j,1)-Image(i,j,2)-Image(i,j,3))/(2*((Image(i,j,1)-Image(i,j,2))^2+(Image(i,j,1)-Image(i,j,3))*(Image(i,j,2)-Image(i,j,3)))^0.5))/pi; end %求H色度 if(Image(i,j,2)>=Image(i,j,3)) HSI(i,j,1) = JD; else HSI(i,j,1) = 360 - JD; end %求S饱和度 HSI(i,j,2) = 1-3*min(min(Image(i,j,1),Image(i,j,2)),Image(i,j,3))/(Image(i,j,1)+Image(i,j,2)+Image(i,j,3)); %求I亮度 HSI(i,j,3) = (Image(i,j,1)+Image(i,j,2)+Image(i,j,3))/3; end end %画图 figure(1); subplot(1,2,1); imshow(uint8(Image)); title('RGB图像') subplot(1,2,2); imshow(uint8(HSI)); title('HSI图像'); end 3.9 HSI to RGB %函数名称为Image_HSI2RGB,输入参数Image,输出参数HSI function Image = Image_HSI2RGB(HSI) %HSI->RGB变换 [lines,samples] = size(HSI(:,:,1)); Image = zeros(lines,samples,3); for i = 1:lines for j = 1:samples if(HSI(i,j,1)>=0 && HSI(i,j,1)=120 && HSI(i,j,1)=240 && HSI(i,j,1)L)变换 JZ = zeros(bands,lines*samples); for i = 1:bands JZ(i,:) = reshape(Image(:,:,i),[1,lines*samples]); end COV = cov(JZ'); [TZXL,TZZ] = eigs(COV); %得到特征值和特征向量 % [TZZ,index] = sort(diag(TZZ),'descend'); %将特征值按降序排列 % TZXL = TZXL(:,index); %特征值为bands列，特征向量矩阵的每一列代表对应特征值的特征向量 A = TZXL'; Y = A * JZ; %Y的每一行代表第n各主成分 end 4 影像融合处理 4.1 基于Brovey变换的影像融合 %函数名称为Image_BroveyChange,输入参数Image1,Image2,输出参数Image3 function Image3 = Image_BroveyChange(Image1,Image2); %Image1表示多光谱图像 [lines1,samples1,bands1] = size(Image1); %Image2表示全色影像 [lines2,samples2,bands2] = size(Image2); Image3 = zeros(lines2,samples2,bands1); for k = 1:bands1 for i = 1:lines2 for j = 1:samples2 Image3(i,j,k) = Image2(i,j,k)*Image1(i,j,k)/sum(Image1(i,j,:)); end end end figure(1) subplot(1,3,1); imshow(uint8(Image1)); title('多光谱影像'); subplot(1,3,2); imshow(uint8(Image2)); title('全色影像'); subplot(1,3,3); imshow(uint8(Image3)); title('Brovey变换融合后的影像'); end 4.2 基于乘积变换的影像融合 %函数名称为Image_MultiplicativeChange,输入参数Image1,Image2,输出参数Image3 function Image3 = Image_MultiplicativeChange(Image1,Image2); %Image1表示多光谱图像 [lines1,samples1,bands1] = size(Image1); %Image2表示全色影像 [lines2,samples2,bands2] = size(Image2); Image3 = zeros(lines2,samples2,bands1); for k = 1:bands1 for i = 1:lines2 for j = 1:samples2 Image3(i,j,k) = Image2(i,j,k)*Image1(i,j,k); end end end Image3 = round((Image3 - min(Image3))*255./(max(Image3) - min(Image3))); %画图 figure(1) subplot(1,3,1); imshow(uint8(Image1)); title('多光谱影像'); subplot(1,3,2); imshow(uint8(Image2)); title('全色影像'); subplot(1,3,3); imshow(uint8(Image3)); title('乘积变换融合后的影像'); end 4.3基于PCA的影像融合 %函数名称为Image_PCAFusion,输入参数Image1,Image2,输出参数Image3 function Image3 = Image_PCAFusion(Image1,Image2) %Image1表示多光谱图像 [lines1,samples1,bands1] = size(Image1); %Image2表示全色影像 [lines2,samples2,bands2] = size(Image2); Image3 = zeros(lines2,samples2,bands1); %PCA影像融合算法 AY = PCA(Image1); %得到PCA系数 A = cell2mat(AY(1)); %得到第一主成分 Y = round(cell2mat(AY(2))); %把高分辨率影像与第一主成分图像做直方图匹配 YY = matching(Image2(:,:,1),reshape(Y(1,:),[lines1,samples1])); %将第一主成分转换成一行 Y(1,:) = reshape(YY,[1,lines2*samples2]); %PCA逆变换 X = double(uint8(inv(A)* Y)); %将矩阵转换成三维矩阵 for i = 1:bands1 Image3(:,:,i) = reshape(X(i,:),[lines1,samples1]); end %画图 figure(1); subplot(1,3,1); imshow(uint8(Image1)); title('高光谱影像'); subplot(1,3,2); imshow(uint8(Image2)); title('全色影像'); subplot(1,3,3); imshow(uint8(Image3)); title('PCA融合后的影像'); end 4.4 基于HSI变换的影像融合 %函数名称为Image_HSIChange,输入参数Image1,Image2,输出参数Image3 function Image = Image_HSIChange(Image1,Image2) %把多光谱图像转换为HSI HSI = RGB2HSI(Image1); %提取出I向量 HSI(:,:,3) = round(HSI(:,:,3)); %把全色影像与I向量做直方图匹配 I = matching(Image2(:,:,1),HSI(:,:,3)); % I = Image2(:,:,1); %把I向量用全色波段替换 HSI(:,:,3) = I; %将替换后的HSI转换到RGB空间 Image3 =HSI2RGB(HSI); figure(1); subplot(1,3,1); imshow(uint8(Image1)); title('高光谱影像'); subplot(1,3,2); imshow(uint8(Image2)); title('全色影像'); subplot(1,3,3); imshow(uint8(Image3)); title('HSI变换融合后的影像'); end 5 影像特征提取 5.1 迭代阈值影像分割算法 function G1 = Image_IterativeThresholdSegmentation(Image) %迭代阈值影像分割算法 Image = Image(:,:,1); [lines,samples] = size(Image); T = median( reshape(Image,[1,lines*samples])); T0 = inf; while( abs(T0 - T) > 1e-2) T0 = T; G1 = Image>=T; K1 = G1.*Image; G2 = Image1,index}(1); y = point{1,index}(2); for i = -1:1 for j = -1:1 xx = x + i; yy = y + j; if( xx>0 && xx0 && yy |
| CopyRight 2018-2019 实验室设备网 版权所有 |