calcaulate_ccsm.pro

idl环境下运行

pro calcaulate_ccsm

rfilename=dialog_pickfile(/read,filter='*.*',title='请选择参考数据')
envi_open_file, rfilename, r_fid=cfid
envi_file_query, cfid, ns=cns, nl=cnl, nb=cnb
cfilename=dialog_pickfile(/read,filter='*.*',title='请选择参考数据')
envi_open_file, cfilename, r_fid=rfid
envi_file_query, rfid, ns=rns, nl=rnl, nb=rnb
map_info = envi_get_map_info(fid=cfid)
ocheck=indgen(cns,cnl,cnb)
cdims = [-1, 0, cns-1, 0, cnl-1]
orefer=indgen(rns,rnl,rnb)
mrefer=indgen(rns,rnl,rnb)
rdims = [-1, 0, rns-1, 0, rnl-1]
pos_arr =indgen(cnb)
;读取refer影像，不移动
  for i=0,cnb-1 do begin
  orefer[*,*,i] = ENVI_GET_DATA(fid=rfid, dims=rdims, pos=pos_arr[i]);获取多波段图像数据到cnb维数组image
  endfor
;  refer=TRANSPOSE(orefer,[2, 0, 1])
  refer=TRANSPOSE(orefer,[2, 0, 1]) ;取消tm反射率的拉伸


;samples=ns ;Hyperion列数
;lines=nl ;Hyperion行数
;bands=nb ;Hyperion波段数


  store=FLTARR(cns,cnl)
  rms_real=FLTARR(cns,cnl,cnb)
  rms_refer=FLTARR(cns,cnl,cnb)
;进行逐个位置的相关系数计算,p代表移动位置
for p=0l,cnb-1 do begin
;   p=2

    ;转换文件读取的起始点位置
    pos_arr=indgen(cnb)
    pos_arr=pos_arr-p
    for i=0,cnb-1 do begin
       if pos_arr[i] lt 0 then  pos_arr[i]=pos_arr[i]+cnb
    endfor
   print,pos_arr
    ;读取文件
    for i=0,cnb-1 do begin
    ocheck[*,*,i] = ENVI_GET_DATA(fid=cfid, dims=cdims, pos=pos_arr[i]);获取多波段图像数据到cnb维数组image
    mrefer[*,*,i] = ENVI_GET_DATA(fid=rfid, dims=rdims, pos=pos_arr[i]);获取多波段图像数据到cnb维数组image
    endfor


    check=TRANSPOSE(ocheck,[2, 0, 1])
    move=TRANSPOSE(mrefer,[2, 0, 1])


   ccsm_real,check=check,refer=refer,rms_real=rms_real,cns=cns,cnl=cnl,cnb=cnb,p=p
   ccsm_refer,check=move,refer=refer,rms_refer=rms_refer,cns=cns,cnl=cnl,cnb=cnb,p=p



endfor

   ENVI_WRITE_ENVI_FILE,rms_real, out_name='rms_real', $
   map_info=map_info

   ENVI_WRITE_ENVI_FILE, rms_refer, out_name='rms_refer', $
   map_info=map_info

  rms_real=TRANSPOSE(rms_real,[2, 0, 1])
  rms_refer=TRANSPOSE(rms_refer,[2, 0, 1])


;       计算交叉相关系数
    for n=0l,cns*cnl-1 do begin
;
     rms_sum=0l

        FOR i =n*cnb,(n+1)*cnb-1 DO begin

        rms_sum=rms_sum+(rms_real[i]-rms_refer[i])^2


        endfor

     store[n]=sqrt(rms_sum/cnb)


     endfor


   ENVI_WRITE_ENVI_FILE, store, out_name='rms_final', $
   map_info=map_info


end