pro statistic
;打开第一幅图像
Ori_Image1=dialog_pickfile(/read,filter='*.*',title='请选择输入基础图像')
envi_open_file,Ori_Image1,r_fid=fid
;查询图像信息
envi_file_query,fid,ns=ns,nl=nl,nb=nb
;将原始图像转为数组，6个波段
arrOri_Image=findgen(ns,nl,2)
dims=[-1,0,ns-1,0,nl-1]
;获得原始数组图像第一个波段数据
arrOri_Image(*,*,0)=envi_get_data(fid=r_fid,dims=dims,pos=0)
;打开第二幅图像
Ori_Image2=dialog_pickfile(/read,filter='*.*',title='请选择第二幅影像')
envi_open_file,Ori_Image2,r_fid=fid
envi_file_query,fid,ns=ns,nl=nl,nb=nb
;获得原始数组图像第二个波段数据
arrOri_Image(*,*,1)=envi_get_data(fid=r_fid,dims=dims,pos=0)
xmean=mean(arrOri_Image(*,*,0))
ymean=mean(arrOri_Image(*,*,1))
result=0D
r_temp2=0D
r_temp4=0D
r_temp6=0D
DI_temp2=0D
g1_temp4=0D
g2_temp4=0D
s1_temp2=0D
s2_temp2=0D
for i=0,ns-1 do begin
for j=0,nl-1 do begin
;计算偏差指数
result1=abs((arrOri_Image(i,j,1)-arrOri_Image(i,j,0))/(arrOri_Image(i,j,0)+0.000000001))
; if arrOri_Image(i,j,0) eq 0 then begin
; result1=0
; endif
result=result+result1
;计算相关系数
r_temp1=(arrOri_Image(i,j,0)-xmean)*(arrOri_Image(i,j,1)-ymean)
r_temp2=r_temp2+r_temp1
r_temp3=(arrOri_Image(i,j,0)-xmean)*(arrOri_Image(i,j,0)-xmean)
r_temp4=r_temp4+r_temp3
r_temp5=(arrOri_Image(i,j,1)-ymean)*(arrOri_Image(i,j,1)-ymean)
r_temp6=r_temp6+r_temp5
;计算光谱扭曲程度
DI_temp1=abs(arrOri_Image(i,j,1)-arrOri_Image(i,j,0))
DI_temp2=DI_temp2+DI_temp1
;计算平均梯度
;第一行第一列差分为0
if i eq 0 or j eq 0 then begin
g1_temp1=0D
g1_temp2=0D
g2_temp1=0D
g2_temp2=0D
endif else begin
g1_temp1=(arrOri_Image(i,j,0)-arrOri_Image(i-1,j,0))*(arrOri_Image(i,j,0)-arrOri_Image(i-1,j,0))
g1_temp2=(arrOri_Image(i,j,0)-arrOri_Image(i,j-1,0))*(arrOri_Image(i,j,0)-arrOri_Image(i,j-1,0))
g2_temp1=(arrOri_Image(i,j,1)-arrOri_Image(i-1,j,1))*(arrOri_Image(i,j,1)-arrOri_Image(i-1,j,1))
g2_temp2=(arrOri_Image(i,j,0)-arrOri_Image(i,j-1,0))*(arrOri_Image(i,j,0)-arrOri_Image(i,j-1,0))
endelse
g1_temp3=sqrt((g1_temp1+g1_temp2)/2)
g1_temp4=g1_temp4+g1_temp3
g2_temp3=sqrt((g2_temp1+g2_temp2)/2)
g2_temp4=g2_temp4+g2_temp3
;计算标准差
s1_temp1=(arrOri_Image(i,j,0)-xmean)*(arrOri_Image(i,j,0)-xmean)
s1_temp2=s1_temp2+s1_temp1
s2_temp1=(arrOri_Image(i,j,1)-ymean)*(arrOri_Image(i,j,1)-ymean)
s2_temp2=s2_temp2+s2_temp1
endfor
endfor
DNmax1=max(arrOri_Image(*,*,0))
DNmin1=min(arrOri_Image(*,*,0))
p1_temp2=make_array(DNmax1)
DNmax2=max(arrOri_Image(*,*,1))
DNmin2=min(arrOri_Image(*,*,1))
p2_temp2=make_array(DNmax2)
p1=0
p2=0
;计算信息熵
for s=DNmin1,DNmax1 do begin
index=where(arrOri_Image(*,*,0)eq s,count1)
p1_temp1=count1/float(ns*nl)
p1_temp2=p1_temp1*(alog10(p1_temp1+0.000000001)/alog10(2))
p1=p1+p1_temp2
endfor
for k=DNmin2,DNmax2 do begin
index=where(arrOri_Image(*,*,1)eq k,count2)
p2_temp1=count2/float(ns*nl)
p2_temp2=p2_temp1*(alog10(p2_temp1+0.000000001)/alog10(2))
p2=p2+p2_temp2
endfor
deviation=result/float(ns*nl)
r=r_temp2/sqrt(r_temp4*r_temp6)
DI=DI_temp2/float(ns*nl)
g1=g1_temp4/float(ns*nl)
g2=g2_temp4/float(ns*nl)
s1=sqrt(s1_temp2/(ns*nl-1))
s2=sqrt(s2_temp2/(ns*nl-1))
print,'第一幅影像均值是',xmean
print,'第二幅影像均值是',ymean
print,'第一幅图像标准差是',s1
print,'第二幅图像标准差是',s2
print,'偏差指数是',deviation
print,'第一幅图像信息熵是',-p1
print,'第二幅图像信息熵是',-p2
print,'第一幅图像平均梯度是',g1
print,'第二幅图像平均梯度是',g2
print,'相关系数是',r
print,'光谱扭曲程度是',DI
end
相关分类