遥感数字图像处理实习报告含Matlab处理代码

更新时间:2023-12-24 15:58:01 阅读量: 教育文库 文档下载

说明:文章内容仅供预览,部分内容可能不全。下载后的文档,内容与下面显示的完全一致。下载之前请确认下面内容是否您想要的,是否完整无缺。

辽宁工程技术大学

《数字图像处理》上机实习报告

教学单位 辽宁工程技术大学 专 业 摄影测量与遥感 实习名称 遥感数字图像处理 班 级 测绘研11-3班 学生姓名 路聚峰 学 号 471120212 指导教师 孙华生

实习1 读取BIP 、BIL、 BSQ文件

一、实验目的

用Matlab读取BIP 、BIL、 BSQ文件,并将结果显示出来。

遥感图像包括多个波段,有多种存储格式,但基本的通用格式有3种,即BSQ、BIL和BIP格式。通过这三种格式,遥感图像处理系统可以对不同传感器获取的图像数据进行转换。BSQ是像素按波段顺序依次排列的数据格式。BIL格式中,像素先以行为单位块,在每个块内,按照波段顺序排列像素。BIP格式中,以像素为核心,像素的各个波段数据保存在一起,打破了像素空间位置的连续性。

用Matlab 读取各个格式的遥感数据,是图像处理的前提条件,只有将图像读入Matlab工作空间,才能进行后续的图像处理工作。

二、算法描述

1.调用fopen函数用指定的方式打开文件。

2.在for循环中调用fread函数,用指定的格式读取各个像素。 3.用reshape函数,重置图像的行数列数。

4.用imadjust函数调整像素的范围,使其有一定对比度。 5.用imshow显示读取的图像。

三、Matlab源代码

1.读取BSQ的源代码: clear all clc lines=400; samples=640; N=6;

img=fopen('D:\\sample_BSQ','rb'); for i=1:N

bi=fread(img,lines*samples,'uint8'); band_cov=reshape(bi,samples,lines);

band_cov2=band_cov'; band_uint8=uint8(band_cov2); tif=imadjust(band_uint8); mkdir('D:\\MATLAB','tifbands1')

name=['D:\\MATLAB\\tifbands1\\tif',int2str(i),'.tif'];

imwrite(tif,name,'tif');

tilt=['波段',int2str(i)];

subplot(3,2,i),imshow(tif);title(tilt); end

fclose(img);

2.读取BIP源代码

clear all

clc lines=400; samples=640; N=6; for i=1:N

img=fopen('D:\\MATLAB\\sample_BIP','rb'); b0=fread(img,i-1,'uint8');

b=fread(img,lines*samples,'uint8',(N-1));

band_cov=reshape(b,samples,lines); band_cov2=band_cov';%×a?? band_uint8=uint8(band_cov2); tif=imadjust(band_uint8); mkdir('E:\\MATLAB','tifbands')

name=['E:\\MATLAB\\tifbands\\tif',int2str(i),'.tif']; imwrite(tif,name,'tif'); %imwrite(A,filename,fmt) tilt=['波段',int2str(i)];

subplot(3,2,i),imshow(tif);title(tilt); fclose(img); end

3.读取BIL的源代码

clear all clc lines=400; samples=640; N=6; for i=1:N

bi=zeros(lines,samples); for j=1:samples

img=fopen('D:\\MATLAB\\sample_BIL','rb');

bb=fread(img,(i-1)*640,'uint8');

b0=fread(img,1*(j-1),'uint8');

bandi_linej=fread(img,lines,'uint8',1*(N*samples-1)); fclose(img);

bi(:,j)=bandi_linej; end

band_uint8=uint8(bi); tif=imadjust(band_uint8); mkdir('D:\\MATLAB','tifbands')

name=['D:\\MATLAB\\tifbands\\tif',int2str(i),'.tif']; imwrite(tif,name,'tif'); tilt=['2¨??',int2str(i)];

subplot(3,2,i),imshow(tif);title(tilt); end

四、运行结果

图1:读取文件的六个波段图

实习2 均值/中值滤波、边缘信息提取

一、实验目的与原理

各种图像滤波算子可以实现图像的增强,去噪,边缘提取等。图像增强的目的在于:1.采用一系列技术改善图像的视觉效果,提高图像的清晰度,2.将图像转换成一种更适合于人或机器进行分析处理的形式。它不是以图像保真度为原则,而是通过处理,设法有选择地突出便于人或机器分析某些感兴趣的信息,抑制一些无用的信息,以提高图像的使用价值。图像增强方法从增强的作用域出发,可分为空间域增强和频率域增强。空间域增强就是直接对图像像素灰度进行操作;频率域增强是对图像经傅里叶变换后的频谱成分进行操作,然后经傅里叶逆变换获得所需结果。

图像滤波可分为空间域滤波和频率域滤波,前者通过窗口或卷积核进行,它参照相邻像素改变单个像素的灰度值。后者对图像进行傅立叶变换,然后对频谱进行滤波。空间域图像滤波称为平滑和锐化,强调像素与其周围相邻像素的关系。去噪滤波为平滑滤波包括均值滤波和中值滤波。锐化滤波包括罗伯特梯度、索伯尔梯度、拉普拉斯算法、定向检测,用以提取线状地物和边缘。此实验用Matlab采用各种滤波对图像进行了处理,处理结果如下: 二、算法描述

1.用imread读取图像文件,并用size获取图像的大小。 2.设计各种滤波算子。

3利用卷积公式对图像的没一个像素进行处理,得到滤波后的图像。 4.用imshow显示滤波后的图像。 三、Matlab源代码 1.均值滤波源码: clear all clc

img=imread('2.jpg');

[row,column,band]=size(img); img0=double(img);

f11=1/9; f12=1/9; f13=1/9; f21=1/9; f22=1/9; f23=1/9; f31=1/9; f32=1/9; f33=1/9;

img1=[img0(:,1,:), img0(:,:,:), img0(:,column,:)]; img2=[img1(1,:,:); img1(:,:,:); img1(row,:,:)]; filtered=zeros(row,column,band); for ii=1: row for jj=1: column

filtered(ii,jj,:)=f11*img2(ii,jj,:) + f12*img2(ii,jj+1,:) + f13*img2(ii,jj+2,:)+ ...

f21*img2(ii+1,jj,:) + f22*img2(ii+1,jj+1,:) + f23*img2(ii+1,jj+2,:) + ...

f31*img2(ii+2,jj,:) + f32*img2(ii+2,jj+1,:) + f33*img2(ii+2,jj+2,:); end end

filtered1=uint8(filtered);

subplot(1,2,1),imshow(img);title('图1 原始RGB图像');

subplot(1,2,2),imshow(filtered1);title('图2 均值滤波后的图像'); imwrite(filtered1,'flower_filtered_mean.jpg'); 2.边缘提取滤波源代码 clear all

img=imread('2.jpg'); [row,column,band]=size(img); img0=double(img); f11=1; f12=0; f13=-1; f21=1; f22=0; f23=-1; f31=1; f32=0; f33=-1;

img1=[img0(:,1,:), img0(:,:,:), img0(:,column,:)]; img2=[img1(1,:,:); img1(:,:,:); img1(row,:,:)]; filtered=zeros(row,column,band); for ii=1: row for jj=1: column

filtered(ii,jj,:)=f11*img2(ii,jj,:) + f12*img2(ii,jj+1,:) + f13*img2(ii,jj+2,:)+ ...

f21*img2(ii+1,jj,:) + f22*img2(ii+1,jj+1,:) + f23*img2(ii+1,jj+2,:) + ...

f31*img2(ii+2,jj,:) + f32*img2(ii+2,jj+1,:) + f33*img2(ii+2,jj+2,:); end end

filtered1=uint8(filtered);

subplot(1,2,1),imshow(img);title('图1 RGB原图像');

subplot(1,2,2),imshow(filtered1);title('图2 边缘提取后的图像'); imwrite(filtered1,'flower_filtered_edge.jpg');

四、运行结果

图1:原始RGB图像 图2:均值滤波后的图像

图3:边缘提取后的图像

实习3 傅里叶变换、傅里叶逆变换,及频域滤波

一、实验目的

按照信号处理理论,根据滤除的频率特征,滤波有3种:1.低通滤波。低通滤波是对频率域的图像通过滤波器H(u,v)削弱或抑制高频部分而保留低频部分的滤波方法。由于图像上的噪声主要集中在高频部分,所以低通滤波可以起到压抑噪声的作用。同时,由于强调了低频成分,图像会变得比较平滑。2.高通滤波。高通滤波是对频率域的图像通过滤波器来突出图像的边缘和轮廓,进行图像锐化

的方法。3.带通滤波。仅保留指定频率范围的滤波,范围外的频率被阻止。

将空间域中的图像变换到频率域中进行计算。空间增强技术强调像元位置和像元之间的关系,但随着考虑的像元数目增多,计算的复杂度增加而且非常耗费计算运算时间,特别是当模板越来越大时,这种现象尤为明显。 频率域增强方法:

1.频率域平滑:保留图像的低频部分而抑制高频部分。 2.频率域锐化:保留图像的高频部分而削弱低频部分。 首先将空间域图像f(x,y)通过傅立叶变换为频率域图像F(u,v),然后选择合适的滤波器H(u,v)对F(u,v)的频谱成分进行增强得到图像G(u,v),再经过傅立叶逆变换将

G(u,v)变换到空间域,得到增强后的图像g(x,y)。

根据傅里叶变换的原理,用Matlab实现对图像的傅里叶变换,再设计各种频率滤波器,包括理想滤波器、巴特沃斯滤波器、指数滤波器等高通或低通滤波器,用这些滤波器对频率图像进行滤波,将得到的滤波后的频率图像经过傅里叶逆变换后得到想要的图像。本实验,对这些滤波器都进行了设计,处理结果如下图: 二、算法描述

1.读取RGB图像,并获得其某个波段。

2.调用fft2对某波段进行傅里叶变换,用fftshift频移函数,使得图像的低频部分移动到频谱的中心。

3.设置低通高通滤波器,对频谱图像进行滤波处理,去除其高频或低频部分。 4.用ifft2对频谱图像进行逆傅里叶变换,得到滤波后的图像。 5.对图像进行归一化,使像素值在0-255之间。 6. imshow显示频谱图像和滤波后的图像。 三、Matlab源代码 1.理想低通滤波源代码 clear all

a=imread('2.jpg'); [X,Y,Z]=size(a);

a1=a(:,:,1); b1= fft2(a1); b2= fftshift(b1); F=abs(b2);

h=zeros(X,Y);

cutoff=0.7; threshold=1-cutoff; lowx=round(X/2-threshold*X/2); upx=round(X/2+threshold*X/2); lowy=round(Y/2-threshold*Y/2); upy=round(Y/2+threshold*Y/2); h(lowx:upx,lowy:upy)=1; lowpass=b2.*h; d0=ifft2(lowpass); result=abs(d0);

result=uint8(result);

F1=log10(reshape(F, X*Y,1));

F2=uint8( ((F1 - min(F1))/(max(F1) - min(F1)))*255); F3=reshape(F2, X, Y);

subplot(1,3,1), imshow(F3), title('Fig.1 ?傅里叶频谱'); subplot(1,3,2), imshow(h), title('Fig.2 理想低通滤波器'); subplot(1,3,3), imshow(result),title('Fig.3 ?理想低通滤波结果');

2.巴特沃斯低通滤波代码: clear all

a=imread('2.jpg'); [X,Y,Z]=size(a);

a1=a(:,:,1);b1= fft2(a1); b2= fftshift(b1); F=abs(b2); h=zeros(X,Y); lowpass=zeros(X,Y); n1=round(X/2); n2=round(Y/2);

dmax=sqrt(n1^2 + n2^2);

cutoff=0.7;d0=(1-cutoff) * dmax; n=3; for x=1:X for y=1:Y

d=sqrt((x-n1)^2+(y-n2)^2); h(x,y)=1/(1+(d/d0)^(2*n)); lowpass(x,y)=b2(x,y).*h(x,y); end end

lp=ifft2(lowpass); result=abs(lp); result=uint8(result); F1=log10(reshape(F, X*Y,1));

F2=uint8( ((F1 - min(F1))/(max(F1) - min(F1)))*255); F3=reshape(F2, X, Y);

subplot(1,3,1), imshow(F3), title('Fig.1 ?傅里叶频谱'); subplot(1,3,2), imshow(h), title('Fig.2 巴特沃斯低通滤波器'); subplot(1,3,3), imshow(result),title('Fig.3 巴特沃斯低通滤波结果');

3.指数低通滤波的主要代码:

for x=1:X for y=1:Y

d=sqrt((x-n1)^2+(y-n2)^2);

h(x,y)=1/(1+exp(-(d/d0))^(-n));

lowpass(x,y)=b2(x,y).*h(x,y);

end end

4.指数高通滤波的主要代码:

for x=1:X for y=1:Y

d=sqrt((x-n1)^2+(y-n2)^2); h(x,y)=1/(1+exp(-(d/d0))^(n)); lowpass(x,y)=b2(x,y).*h(x,y); end end

四、运行结果

图1:傅里叶频谱

图2:理想低通滤波器 图3:傅里叶低通滤波结果

图4:理想高通滤波器 图5:傅里叶高通滤波结果

图6:巴特沃斯低通滤波器 图7:巴特沃斯低通滤波结果

图8:巴特沃斯高通滤波器 图9:巴特沃斯高通滤波结果

图10:指数高通滤波器 图11:指数高通滤波结果

图12:指数低通滤波器 图13:指数低通滤波结果

实习4 主成分变换、主成分逆变换

一、实验目的与原理

遥感多光谱影像波段多,信息量大,对图像解译很有价值。但数据量太大,在图像处理计算时,也常常耗费大量的机时,占据大量的磁盘空间。实际上,一些波段的遥感数据之间都有不同程度的相关性,存在着数据冗余。多光谱变换方法可通过函数变换,达到保留主要信息、降低数据量、增强或提取有用信息的目的。其变换的本质是对遥感图像进行线性变换,使多光谱空间的坐标系按一定规律进行旋转。多光谱变换主要有两种变换:主成分变换和缨帽变换。 (1)主成分变换(K-L变换)的目的:

1.K-L变换是图像分析与模式识别中的重要工具,用于特征抽取,降低特2.征数据的维数

3.K-L变换用于图像压缩时可以实现有损压缩和无损压缩

4.K-L变换后的N幅图像统计上互不相关,因此K-L变换可以去除图像数据的相关性,提取主要信息。

用Matlab实现六个波段图像的主成分分析。求出六个波段的协方差阵,再求协方差阵的特征值、特征向量。只取特征值比较大的数,其占各个波段主要成分,特征值小的,大部分噪声等内容。 二、算法描述

主成分变换过程:

1.读取六个波段的图像。

2.将每个波段的所有像素存成一列,共六列,形成MN*6矩阵。 3.求出每个波段的平均像素值。

4.利用算出的平均像素值,求MN*6矩阵的协方差阵。 5.求协方差阵的特征值与特征向量。

6.将特征值按从大到小排列,特征向量与特征值对应。 主成分的逆变换过程: (1)确定保留的波段数

(2)将去除的波段用0值替代

(3)将去噪后的PCA结果 SCORE*inv( COEFF) (4)去中心化(各波段加上其均值) (5)重新排列图像行列数 三、Matlab源代码 1.主成分分析代码 clear all clc

num_band=6; %% 2¨??êy b1=imread('tif1.tif'); b2=imread('tif2.tif'); b3=imread('tif3.tif'); b4=imread('tif4.tif'); b5=imread('tif5.tif'); b6=imread('tif6.tif'); [x,y]=size(b1);

img0=[reshape(b1,x*y,1),reshape(b2,x*y,1),

reshape(b3,x*y,1),reshape(b4,x*y,1),reshape(b5,x*y,1),reshape(b6,x*y,1)];

img1=double(img0); average=mean(img1);

img3=img1-repmat(average,x*y,1); c=cov(img3);

[coff,lam] = eig(c); lam1=lam(6:-1:1,6:-1:1); coff2=coff(:,6:-1:1); lam_val=diag(lam1); lam_sum=sum(lam_val); lam_n_sum=sum(lam_val(1:3)); prop=lam_n_sum/lam_sum; result=img3*coff2;

result2=reshape(result,x,y,6);

result3=uint8(result2(:,:,1)-min(min(result2(:,:,1)))/max(max(result2(:,:,1))-min(min(result2(:,:,1))))*255); imshow(result3); hold on

result3=uint8(result2(:,:,2)-min(min(result2(:,:,2)))/max(max(result2(:,:,2))-min(min(result2(:,:,2))))*255); figure,imshow(result3);

result3=uint8(result2(:,:,3)-min(min(result2(:,:,3)))/max(max(result2(:,:,3))-min(min(result2(:,:,3))))*255); figure,imshow(result3);

result3=uint8(result2(:,:,4)-min(min(result2(:,:,4)))/max(max(result2(:,:,4))-min(min(result2(:,:,4))))*255); figure,imshow(result3);

result3=uint8(result2(:,:,5)-min(min(result2(:,:,5)))/max(max(result2(:,:,5))-min(min(result2(:,:,5))))*255); figure,imshow(result3);

result3=uint8(result2(:,:,6)-min(min(result2(:,:,6)))/max(max(result2(:,:,6))-min(min(result2(:,:,6))))*255); figure,imshow(result3); result2(:,:,4:6)=0;

result_new=reshape(result2,x*y,6); result_ipca=result_new*inv(coff2);

result_ipcaz2=result_ipca+repmat(average,x*y,1); result_ipcaz3=reshape(result_ipcaz2,x,y,6); figure,imshow(uint8(result_ipcaz3(:,:,1))); C = cat(3, uint8(result_ipcaz3(:,:,3)),

uint8(result_ipcaz3(:,:,4)),uint8(result_ipcaz3(:,:,5))); figure,imshow(C);

四、运行结果

图1:第一主分量 图2:第二主分量

图3 :第三主分量 图4:第四主分量

图5:第五主分量 图6:第六主分量

图7:逆变换后的图像 图8:逆变换后的彩色合成图像

由图可以看出,图像的主要成分主要集中在前三个主分量图像上,后三个图像大部分都是噪声。通过逆变换恢复原来的图像,这样,经过K-L变换后的前三个主分量,信噪比大,噪声相对小,突出了主要信息,达到了图像增强的目的。另外,经过K-L变换后,第一或前二或前三个主分量已经包含了绝大多数的地物信息,足够分析使用,同时数据量大大减少。应用中,常常只取前三个主分量作假彩色合成(如图8),数据量可减少到43%,既实现了数据压缩,也可作为分类前的特征选择。

实习5 缨帽变换

一、 实验目的

缨帽变换(K-T变换)的目的:该变换也是一种坐标空间发生旋转的线性组合变换,但旋转后的坐标轴不是指向主成分方向,而是指向与地物特别是和植被生成以及土壤有密切关系的方向 二、 算法描述

1. 读取六个波段的图像。

2. 将每个波段的所有像素排成一列,共六列组成一个矩阵。 3. 设计转换矩阵。

4. 图像矩阵与转换矩阵相乘,得到六个相互垂直的分量。 5. 逐一显示各个分量。 三、Matlab源代码 clear all clc

b1=imread('tif1.tif'); b2=imread('tif2.tif'); b3=imread('tif3.tif'); b4=imread('tif4.tif'); b5=imread('tif5.tif'); b6=imread('tif6.tif'); [x,y]=size(b1);

img0=[reshape(b1,x*y,1),reshape(b2,x*y,1),

reshape(b3,x*y,1),reshape(b4,x*y,1),reshape(b5,x*y,1),reshape(b6,x*y,1)];

img0=double(img0);

KT_tm = [0.3561, 0.3972, 0.3904, 0.6966, 0.2286, 0.1596; -0.3344,-0.3544,-0.4556, 0.6966,-0.0242,-0.2630; 0.2626, 0.2141, 0.0926, 0.0656,-0.7629,-0.5388; 0.0805,-0.0498, 0.1950,-0.1327, 0.5752,-0.7775; -0.7252,-0.0202, 0.6683, 0.0631,-0.1494,-0.0274; 0.4000,-0.8172, 0.3832, 0.0602,-0.1095, 0.0985]; img1 = img0 * KT_tm'; min0 = min(img1); max0 = max(img1); dm = max0 - min0;

img11 = img1 - repmat(min0,x*y,1); img2 = (img11./repmat(dm,x*y,1))*255; image = reshape(img2,x,y,6);

figure,imshow(uint8(image(:,:,1))),title('KT变换第一分量:亮度'); figure,imshow(uint8(image(:,:,2))),title('KT变换的第二分量:绿度'); figure,imshow(uint8(image(:,:,3))),title('KT变换的第三分量:湿度'); figure

imshow(uint8(image(:,:,4))),title('KT变换的第四分量:黄度指数及噪声');

四、运行结果

图9:第一分量:亮度图 图10:第二分量:绿度

图11:第三分量:湿度 图12:第四分量:黄度指数及噪声 用Matlab编写代码实现缨帽变换,其基本思想就是:采用变换矩阵,使图像指向与地物特别是和植被生长以及土壤有密切关系的方向。变换后

Y?[y1,y2,y3,y4,y5,y6],这六个分量相互垂直,前三个分量具有明确的地物意义:y1为亮度,实际上是TM六个波段的加权和,反映出图像总体的反射值;y2为绿度,y3为湿度。

实习6 图像分类

(最小距离分类、最大似然分类、K均值分类)

一、实验目的与原理

图像分类的目的是将图像中每个像元根据其不同波段的光谱亮度、空间结构特征或者其他信息,按照某种规则或算法划分为不同的类别。遥感图像分类就是对地球表面及其环境在遥感图像上的信息进行识别和分类,从而达到识别图像信息所对应的实际地物,提取所需地物信息的目的。

遥感图像分类的原理:由于地物的成分、性质、分布情况的复杂程度和成像条件不同,以及一个像元或瞬时试场里往往有两种或多种地物的情况,即混合像元,使得同类地物的特征向量也不尽相同,而且使得不同地物类型的特征向量之间的差别也不都是显著的。遥感图像的光谱特征通常是以地物在多光谱图像上的灰度体现出来的,即不同地物在同一波段图像上表现的灰度一般互不相同;同时,不同地物在多个波段上的灰度呈现规律也不相同,这就构成了我们在图像上区分不同地物的物理依据。

遥感图像传统的计算机分类方式有监督分类和非监督分类两种。监督分类包括最小距离分类法、平行多面体分类法、最大似然分类法。非监督分类包括K-

均值聚类法、ISODATA分类法。本实验将用Matlab实现监督分类的最小距离分类和最大似然分类,以及非监督分类中的K-均值分类。 三、 算法描述

(1) 最小距离分类算法描述

1.读取各个波段文件,并加载各类地物在各个波段的均值文件。 2.计算每个像元向量到各类地物的距离。

3.比较像元到各类地物的距离,取最小的距离,并将像元归为该类。 4.为各类赋予一定的颜色。 (2)最大似然分类算法描述:

1.读取BIP文件,将每个波段读为一列。 2.加载各类地物在各个波段的均值和方差

3.对每个像元向量计算属于各个类的概率,比较概率的大小,概率最大的将像元归为该类。

4.对每个类赋予一定的颜色。 5.显示分类后的图像。

(3)K-means分类算法描述:

1.为每个聚类确定一个初始的聚类中心(一定要在最小值和最大值之间,这一

步非常重要)

2.将样本集中的每一个样本按照最小距离原则,分配到k个聚类中的某一个 3.使用每个聚类中所有样本的均值作为新的聚类中心

4.如果聚类中心有变化则重复2、3步直到聚类中心不再变化为止 5.最后得到的k个聚类中心就是聚类的结果 四、 Matlab源代码

(1)最小距离分类代码: clear all clc

x = 400; y = 640; N = 6;

image = zeros(x*y,N); for i = 1:N

img=fopen('sample_BIP','rb'); b0 = fread(img,i-1,'uint8');

image(:,i) = fread(img,x*y,'uint8',N-1); fclose(img); end

load crop.txt; load forest.txt; load water.txt; load soil1.txt; load soil2.txt; load soil3.txt; dis = zeros(x*y,N);

dis(:,1) = sqrt(sum((image-repmat(crop(1,:),x*y,1)).^2,2));

dis(:,2) = sqrt(sum((image-repmat(forest(1,:),x*y,1)).^2,2)); dis(:,3) = sqrt(sum((image-repmat(water(1,:),x*y,1)).^2,2)); dis(:,4) = sqrt(sum((image-repmat(soil1(1,:),x*y,1)).^2,2)); dis(:,5) = sqrt(sum((image-repmat(soil2(1,:),x*y,1)).^2,2)); dis(:,6) = sqrt(sum((image-repmat(soil3(1,:),x*y,1)).^2,2)); min_dis = min(dis,[],2); classes = zeros(x*y,3); for i = 1:x*y

ind = find(dis(i,:)==min_dis(i)); switch ind case 1

classes(i,:) = [255,0,0]; case 2

classes(i,:) = [0,255,0]; case 3

classes(i,:) = [0,0,255]; case 4

classes(i,:) = [255,255,0]; case 5

classes(i,:) = [255,0,255]; case 6

classes(i,:) = [0,255,255]; end end

img2 = reshape(classes,y,x,3);

RGB_classes = cat(3,img2(:,:,1)',img2(:,:,2)',img2(:,:,3)'); imshow(uint8(RGB_classes)); title('最小距离分类');

(3)最大似然分类代码: clear all clc

x = 400; y = 640; N = 6;

image = zeros(x*y,N); for i = 1:N

img=fopen('sample_BIP','rb'); b0 = fread(img,i-1,'uint8');

image(:,i) = fread(img,x*y,'uint8',N-1); fclose(img); end

load crop.txt; load forest.txt; load water.txt;

load soil1.txt; load soil2.txt; load soil3.txt;

all=cat(3,crop,forest,water,soil1,soil2,soil3); cov = zeros(6,6,6); mean = zeros(1,6,6); for i = 1:6

mean(1,:,i) = all(1,:,i); cov(:,:,i)=diag(all(2,:,i)); end

dk = zeros(x*y,6); for i = 1:x*y for j = 1:6

dk(i,j) = -log(det(cov(:,:,j)))*0.5 - 0.5*(image(i,:) - mean(1,:,j))/cov(:,:,j) * (image(i,:) - mean(1,:,j))'; end end

dk_max = max(dk,[],2); classes = zeros(x*y,3); for i = 1:x*y

ind = find(dk(i,:)==dk_max(i)); switch ind case 1

classes(i,:) = [255,0,0]; case 2

classes(i,:) = [0,255,0]; case 3

classes(i,:) = [0,0,255]; case 4

classes(i,:) = [255,255,0]; case 5

classes(i,:) = [255,0,255]; case 6

classes(i,:) = [0,255,255]; end end

img2 = reshape(classes,y,x,3);

RGB_classes = cat(3,img2(:,:,1)',img2(:,:,2)',img2(:,:,3)'); imshow(uint8(RGB_classes)); title('最大似然分类'); (3)K-means分类代码: clear all clc x = 400;

y = 640; N = 6;

image = zeros(x*y,N); for i = 1:N

img=fopen('sample_BIP','rb'); b0 = fread(img,i-1,'uint8');

image(:,i) = fread(img,x*y,'uint8',N-1); fclose(img); end

load ini_center.txt;

center = zeros(size(ini_center)); threshold = 1e-3; max_d = 1;

while max_d > threshold dis = zeros(x*y,N);

dis(:,1) = sqrt(sum((image-repmat(ini_center(1,:),x*y,1)).^2,2)); dis(:,2) = sqrt(sum((image-repmat(ini_center(2,:),x*y,1)).^2,2)); dis(:,3) = sqrt(sum((image-repmat(ini_center(3,:),x*y,1)).^2,2)); dis(:,4) = sqrt(sum((image-repmat(ini_center(4,:),x*y,1)).^2,2)); dis(:,5) = sqrt(sum((image-repmat(ini_center(5,:),x*y,1)).^2,2)); dis(:,6) = sqrt(sum((image-repmat(ini_center(6,:),x*y,1)).^2,2)); min_dis = min(dis,[],2); classes = zeros(x*y,3); crop = zeros(1,6); forest = zeros(1,6); water = zeros(1,6); soil1 = zeros(1,6); soil2 = zeros(1,6); soil3 = zeros(1,6); n_crop = 0; n_forest = 0; n_water = 0; n_soil1 = 0; n_soil2 = 0; n_soil3 = 0; for i = 1:x*y

ind = find(dis(i,:) == min_dis(i)); switch ind(1,1) case 1

crop = crop + image(i,:); n_crop = n_crop + 1; classes(i,:) = [255,0,0]; case 2

forest = forest + image(i,:);

n_forest = n_forest + 1; classes(i,:) = [0,255,0]; case 3

water = water + image(i,:); n_water = n_water + 1; classes(i,:) = [0,0,255]; case 4

soil1 = soil1 + image(i,:); n_soil1 = n_soil1 + 1; classes(i,:) = [255,255,0]; case 5

soil2 = soil2 + image(i,:); n_soil2 = n_soil2 + 1; classes(i,:) = [255,0,255]; case 6

soil3 = soil3 + image(i,:); n_soil3 = n_soil3 + 1; classes(i,:) = [0,255,255]; end end

center(1,:) = crop./n_crop; center(2,:) = forest./n_forest; center(3,:) = water./n_water; center(4,:) = soil1./n_soil1; center(5,:) = soil2./n_soil2; center(6,:) = soil3./n_soil3;

max_d = max(max(abs(center-ini_center))); ini_center = center; end

img2 = reshape(classes,y,x,3);

RGB_classes = cat(3,img2(:,:,1)',img2(:,:,2)',img2(:,:,3)'); imshow(uint8(RGB_classes)),title('kmeans·?àà?á1?');

四、运行结果

图1:最小距离分类

图2:最大似然分类

图3:K-means分类

实习7 大气校正、反射率、地表亮温反演

实习8 Habib教授摄影测量课程总结

主要讲述了摄影测量里的共线方程。由共线方程求坐标的正算与反算。即由内方位元素、外方位元素、像点坐标去求解地面点坐标,及由内方位元素、地面点坐标以及相应的像点坐标去求外方位元素。并考虑到影像扭曲变形的影响,采用一定的算法消除。

x?xp?cy?yp?cr11?(X?Xo)?r21?(Y?Yo)?r31?(Z?Zo) (1)

r13?(X?Xo)?r23?(Y?Yo)?r33?(Z?Zo)r12?(X?Xo)?r22?(Y?Yo)?r32?(Z?Zo) (2)

r13?(X?Xo)?r23?(Y?Yo)?r33?(Z?Zo)

图1:共线方程示意图

R(?,?,?)(X0,Y0,Z0)投影中心在大地坐标系中的坐标,与地面点同一坐标系,为旋转矩阵,即外方位元素。a(x,y)为像点坐标。

Habib教授还主要讲述了两个主要的Matlab程序。第一个是由内外方位元素以及像点坐标去求地面点坐标,其中考虑到了影像扭曲变形的影响,并采用一定的算法消除。第二个是由内方位元素,地面控制点坐标,以及近似外方位元素,通过最小二乘迭代法去求外方位元素。 求地方点坐标的代码主要实现步骤:

1. 首先,由外方位元素?,?,?,求解旋转矩阵,定义旋转矩阵函数function [R] =

Rotation(W, P, K)。

2. 由扭曲形变改正算法,定义一个形变函数,用来求解形变量,function [dist_x,

dist_y] = Estimate_Distortion(x, y, xp, yp, k1, k2, k3, p1, p2, p3, A1, A2)。

3. 将内方位元素,形变改正系数,地面点坐标的文件,加载到工作空间。 4. 通过旋转矩阵,将经过形变改正的左右影像的像点坐标转换到地面坐标系

下,使其与地面点坐标在同一坐标系下。得到VL、VR两个向量。

5. 由左右影像的外方位元素(XO_L,YO_L,ZO_L)和(XO_R,YO_R,ZO_R)作

差求出基线向量(BX,BY,BZ)。

?VL(1)6. ??VL(2)??VL(3)?VR(1)??BX?X(1)???BY?,通过解方程,求的左右向量的比例系数

?VR(2)????X(2)???????VR(3)???BZ???X(1)??X(2)?。 ??这样便得地面点坐标:XG_L = XO_L + X(1) * VL(1); YG_L = YO_L + X(1) * VL(2); ZG_L = ZO_L + X(1) * VL(3);

7.

本文来源:https://www.bwwdw.com/article/e9u5.html

微信扫码分享

《遥感数字图像处理实习报告含Matlab处理代码.doc》
将本文的Word文档下载到电脑,方便收藏和打印
推荐度:
点击下载文档
下载全文
范文搜索
下载文档
Top