实验2_纹理_结构特征分析一、实验目的加深对SAR影像纹理、结构特征及相关概念的理解，能编程实现SAR影像纹理、结构特征的统计分析。

二、实验软件Matlab三、实验数据德国TerraSAR卫星获取的成都市部分SAR影像sar_cd01.tif及从中裁剪的两部分影像sar_cd01_1.tif、sar_cd01_2.tif。

四、实验内容及步骤1. 基于Matlab，编写统计分析影像纹理特征和结构特征的程序；1.1编写纹理特征程序1）纹理特征主程序CTFmain：依次实现的功能为：a打开并读取影像文件；b获取灰度影像；c扩展原始数据；d统计纹理特征频次；e对纹理特征频次排序并得到原始编码索引；f要求用户输入纹理特征分类个数目；g对纹理特征分类并得到类别编码对原始编码的索引；h生成纹理特征图。

function im=CTFmain% 选择影像文件[fn,pn]=uigetfile({'*.jpg;*.tif;*.png;*.gif',...'All Image Files(*.jpg,*.tif,*.png,*.gif)'; ...'*.fig','Figures(*.fig)';'*.*','All Files' },...'请选择影像文件');if fn~=0% 添加当前路径addpath(pn);% 读取影像文件dat=imread(fn);% 取出起始时间t1=clock;% 获取灰度影像data0=GetGMap(dat);fprintf('%s\n','灰度影像生成完毕！');% 扩展原始数据矩阵data0为datadata=EnlargeMat(data0,1,1);fprintf('%s\n','影像数据扩展完毕！');% 统计每种纹理特征向量出现的频数，得到纹理特征查询矩阵stvm和原始编码图codm [stvm,codm]=StatTVm(data);fprintf('%s\n','纹理特征统计完毕！');% 对各种纹理特征向量频数由大到小排序，得到排序后的纹理特征原始编码索引id0=SortGetId0(stvm);% 取出终了时间t2=clock;% 显示程序运行时间fprintf('%s%.1f%s\n','该阶段运行时间：',etime(t2,t1),'s');% 提示用户输入纹理特征分类数目n=input('请根据统计数据输入纹理特征类别的数目n [10]:');if isempty(n)n=10;end% 取出起始时间t1=clock;% 合并相似的纹理特征，得到频次和类别编码索引fid以及原始编码和类别编码索引iid [fid,iid]=IcprtTF(stvm,id0,n);fprintf('%s\n','特征类别编码完毕！');% 生成纹理特征图im=CreatIM(codm,fid,iid);fprintf('%s\n','纹理特征图生成完毕！');% 取出终了时间t2=clock;% 显示程序运行时间fprintf('%s%.1f%s\n','该阶段运行时间：',etime(t2,t1),'s');% 显示纹理特征图imshow(im,[min(min(im)),max(max(im))]);elseerror('您没有选择影像文件！');endend2）获取灰度图像子程序GetGMap：实现了根据输入数据结构采用相应算法得到灰度图像的功能，即：a如果是单层影像数据，直接认为是灰度图像，有data0 = dat；b如果是3层影像数据，采用gray = 0.299*r+0.587*g+0.114*b公式计算出灰度图像：c如果是多余3层的影像数据，采用所有数据层平均的方法得到灰度图像。

function data0=GetGMap(dat)ds=size(dat);data0=zeros(ds(1),ds(2));dss=size(ds);if dss(2)==3data0=0.299*dat(:,:,1)+0.587*dat(:,:,2)+0.114*dat(:,:,3);else if dss(2)>3for i=1:dss(2)data0=data0+dat(:,:,i);endelsedata0=dat;endendend3）扩展原始数据子程序EnlargeMat：实现将数据data0扩展为data的功能，具体的方法是：a将data0的4个角点值分别赋给在其所处对角线另一端的data的角点；b将data0的初始端制定数目的行赋给data最末端制定数目的行，而data0最末端制定数目的行赋给data初始端制定数目的行；c与b相同的方法给data的初始端制定数目的列和最末端制定数目的列赋值。

function data=EnlargeMat(data0,dr,dc)% 取出原始数据的行列数[row,col]=size(data0);% 扩展data0为datadata=zeros(row+2*dr,col+2*dc);data(dr+1:row+dr,dc+1:col+dc)=data0;data(1,1)=data0(row,col);data(row+2*dr,1)=data0(1,col);data(1,col+2*dc)=data0(row,1);data(row+2*dr,col+2*dc)=data0(1,1);data(1,dc+1:col+dc)=data0(row,:);data(row+2*dr,dc+1:col+dc)=data0(1,:);data(dr+1:row+dr,1)=data0(:,col);data(dr+1:row+dr,col+2*dc)=data0(:,1);end4）统计纹理特征子程序StatTVm：实现了对每种纹理特征出现频次的统计，将每种纹理特征向量按三进制转换为十进制再加1作为其对应的原始编码，然后删除纹理特征向量频次统计矩阵（stvm）中频次为0的纹理特征所对应的行数据，最后返回stvm和赋予了原始编码的原始编码图（codm）。

function [stvm,codm]=StatTVm(data)% 取出输入数据data的行列数[row,col]=size(data);% 定义纹理特征向量统计矩阵stvmstvm=CreatTVm;% 定义纹理特征编码图codm=zeros(row-2,col-2);% 初始化一个进度条h=waitbar(0,'正在统计纹理特征：请稍后...');% 计算tv0，并对比stvm中每一种纹理特征向量，统计其出现的次数for i=2:row-1for j=2:col-1w=data(i-1:i+1,j-1:j+1);tv=ComputTV(w);k=ComputIdx(tv);for k0=1:6561if k==k0stvm(k,9)=stvm(k,9)+1;codm(i-1,j-1)=k0;endendendwaitbar((i-1)/(row-1));endid=stvm(:,9)==0;stvm(id,:)=[];waitbar((row-1)/(row-1));% 关闭进度条close(h);end5）生成纹理特征向量频次统计矩阵子函数CreatTVm：实现生成一个纹理特征向量频次统计矩阵的功能，其最后两列存频数和原始编码。

function tvm=CreatTVm% 定义纹理特征向量矩阵，最后两列存频数和原始编码tvm=zeros(6561,10);% 生成纹理特征向量矩阵，其包含了6561种纹理特征向量i=1;for k1=0:2for k2=0:2for k3=0:2for k4=0:2for k5=0:2for k6=0:2for k7=0:2for k8=0:2tvm(i,1)=k1;tvm(i,2)=k2;tvm(i,3)=k3;tvm(i,4)=k4;tvm(i,5)=k5;tvm(i,6)=k6;tvm(i,7)=k7;tvm(i,8)=k8;tvm(i,10)=i;i=i+1;endendendendendendendendend6）计算纹理特征向量子函数ComputTV：实现根据3*3窗口w的数据计算对应的纹理特征向量的功能。

function tv=ComputTV(w)% 定义纹理特征向量tvtv=zeros(1,9);% 计算纹理特征向量tvfor r=1:3for c=1:3k=3*r+c-3;if w(r,c)==w(2,2)tv(1,k)=1;else if w(r,c)>w(2,2)tv(1,k)=2;elsetv(1,k)=0;endendendend% 删除w窗口中心所对应的列tv(:,5)=[];end7）计算原始编码子函数ComputIdx：实现根据纹理特征向量计算出其对应的原始编码的功能。

function idx=ComputIdx(tv)% 计算纹理特征向量tv的索引值idx=tv(8)+3*tv(7)+3^2*tv(6)+3^3*tv(5)+3^4*tv(4)+3^5*tv(3)+3^6*tv(2)+3^7*tv(1)+1;end8）建立纹理特征频次排序索引子程序SortGetId0：实现了对精简后的纹理特征频次由大到小排序，得到排序后的纹理特征原始编码索引的功能。

function id0=SortGetId0(stvm)% 取出stvm大小s=size(stvm);% 对频数排序[fqn,id0]=sort(stvm(:,9),'descend');% 给id0替换原始编码值for i=1:s(1)jd=id0(:)==i;id0(jd)=stvm(i,10);endend9）纹理特征合并分类编码子程序IcprtTF：实现了合并相似的纹理特征，得到频次和类别编码索引fid以及原始编码和类别编码索引iid的功能。

function [fid,iid]=IcprtTF(stvm,id0,n)% 取出stvm的大小s=size(stvm);% 初始化频数和类别的索引fidfid=zeros(n,2);% 初始化编码和类别的索引iidiid=zeros(s(1),2);% 初始化频数较大的前n种纹理特征查询表stv0stv0=zeros(n,9);% 给stv0赋值，前8位是纹理特征向量，最后1位是纹理特征对应的编码for k=1:nfor i=1:s(1)if id0(k)==stvm(i,10)stv0(k,1:8)=stvm(i,1:8);stv0(k,9)=id0(k);endendend% 将剩余向量和查询表里的向量按夹角归类,建立编码和类别的索引iid for i=1:s(1)ag=pi*ones(n,1);for k=1:nif stvm(i,10)~=stv0(:,9)ag(k)=ComputAg(stv0(k,1:8),stvm(i,1:8));elsebreak;endend[agm,k0]=min(ag);if agm==piiid(i,1)=stvm(i,10);iid(i,2)=stvm(i,10);elseiid(i,1)=stvm(i,10);iid(i,2)=stv0(k0,9);endend% 给fid赋值for k=1:nfor i=1:s(1)if iid(i)==id0(k)fid(k,1)=fid(k,1)+stvm(i,9);endendendfid(:,2)=id0(1:n);end10）生成纹理特征图子程序CreatIM：实现生成纹理特征图的功能：首先建立灰度查找表glt（Gray Look Table），然后根据原始编码图codm和频次-类别编码索引fid和原始编码-类别编码iid的信息，在glt中查找对应灰度，给纹理特征强度图im（intensity map）赋值。