当前位置:首页 >> 农林牧渔 >>

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


辽宁工程技术大学

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

教学单位 专 业

辽宁工程技术大学 摄影测量与遥感 遥感数字图像处理 测绘研 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';%×??? 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=['?¨??',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,:) + f13*img2(ii,jj+2,:)+ ... f21*img2(ii+1,jj,:) f23*img2(ii+1,jj+2,:) + ... f31*img2(ii+2,jj,:) 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,:) + f13*img2(ii,jj+2,:)+ ... + f12*img2(ii,jj+1,:) + f32*img2(ii+2,jj+1,:) + + f22*img2(ii+1,jj+1,:) + + f12*img2(ii,jj+1,:)

f21*img2(ii+1,jj,:) f23*img2(ii+1,jj+2,:) + ... f31*img2(ii+2,jj,:) f33*img2(ii+2,jj+2,:); end end filtered1=uint8(filtered);

+ +

f22*img2(ii+1,jj+1,:) f32*img2(ii+2,jj+1,:)

+ +

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; %% ?¨???? 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 , y 2 , y 3 , y 4 , y 5 , y 6 ] ,这六个分量相互垂直,前三个分量具有明确的地物意义:

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:最小距离分类

图2:最大似然分类

图3:K-means分类

实习 7

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

实习 8

Habib 教授摄影测量课程总结

主要讲述了摄影测量里的共线方程。由共线方程求坐标的正算与反算。即由 内方位元素、外方位元素、像点坐标去求解地面点坐标,及由内方位元素、地面 点坐标以及相应的像点坐标去求外方位元素。并考虑到影像扭曲变形的影响,采 用一定的算法消除。
x ? xp ?c y ? yp ?c r11 ? ( X ? X o ) ? r 21 ? (Y ? Y o ) ? r 31 ? ( Z ? Z o ) r13 ? ( X ? X o ) ? r 23 ? ( Y ? Y o ) ? r 33 ? ( Z ? Z o ) r12 ? ( X ? X o ) ? r 22 ? (Y ? Y o ) ? r32 ? ( Z ? Z o ) r13 ? ( X ? X o ) ? r 23 ? (Y ? Y o ) ? r33 ? ( Z ? Z o )

(1) (2)

图 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)。 6.
? V L (1) ? VL (2) ? ? V L (3) ? ? V R (1) ? ? ?V R ( 2 ) ? ? V R (3) ? ? ? BX ? X (1) ? ? ? ? ? ? BY ? X (2)? ? BZ ? ? ? ,通过解方程,求的左右向量的比例系数 ? ? ?

? X (1) ? ? ?。 ? X (2)?

这样便得地面点坐标:XG_L = XO_L + X(1) * VL(1); ZG_L = ZO_L + X(1) * VL(3);

YG_L = YO_L + X(1) * VL(2);

7.


相关文章:
遥感数字图像处理教程实习报告
5 《数字图像处理实习报告 实习二: 实习二:灰度变换实习一、实验目的在学习遥感数字图像灰度变换理论基础上,应用所学知识,基于遥感图像处 理软件 ENVI 进行图像...
遥感数字图像处理2
遥感数字图像处理2_天文/地理_自然科学_专业资料。遥感数字图像处理实习报告 实习序号及题目 实习人姓名 专业班级 实习地点 影像几何精纠正实习指导教师姓名 地址 ...
遥感数字图像处理实验报告
遥感数字图像处理实验报告_实习总结_总结/汇报_实用文档。遥感数字图像处理 学院 班级 学号 姓名 理学院 地信 131 编写日期:2015.5 1 ?? 作业 a 1.LS8_C_...
遥感数字图像处理理论课大纲
归纳总结等为主要实验方法,以实验报告、 论文等为最终形式,以学生为主体的实践...ENVI 遥感图像处理及 MATLAB 软件熟悉,2; 2. 遥感数据的预处理,2; 3. 遥感...
遥感数字图像处理3
遥感数字图像处理3_天文/地理_自然科学_专业资料。遥感数字图像处理实习报告 实习序号及题 目 实习人姓名 e-mail 地址 实习地点 遥感影像反差增强专业班级 实习指导...
遥感数字图像处理-matlab-频域处理
遥感数字图像处理-matlab-频域处理_IT/计算机_专业资料...图像填充掌握折叠误差产生的原因及处理方法,以矩形...地理学与遥感科学学院 蓝建航 ①实验代码 { close...
遥感数字图像处理-matlab-主成份及穗帽变换
遥感数字图像处理-matlab-主成份穗帽变换_工学_高等...EVNI raw 格式图像的读取和显示①实验代码 { close...
遥感图像处理实习总结
遥感实习总结 专业:摄影测量与遥感技术 班级: 姓名: 学号: 为期两周的遥感数字图像处理结束了,在老师的精心安排 下,我们全身心的投入到这次实习中。虽然是满天的...
遥感图像处理实验报告
遥感图像处理实验报告_实习总结_总结/汇报_应用文书。使用ENVI软件处理遥感数字图像《遥感数字图像处理实习报告 学院: 环境与资源学院 班级: 学号: 姓名: 地理 ...
遥感数字图像处理教程整理
遥感数字图像处理教程整理_工学_高等教育_教育专区。辐射分辨率:辐射分辨率是指传感器...辐射校正包括 3 部分的内容:传感器端的辐射校正、大气校正和地表辐射校正 ? ...
更多相关标签: