fx4.pro

来自「用分形算法对高光谱图像压缩的MATLAB程序」· PRO 代码 · 共 173 行

pro fx4,xll
t0=systime(1)
;xll=indgen(48,48,48)
DSijt=fltarr(4,4,4)      ;缩小后的D数组  dblarr
cens=intarr(12,12,12)
hans=intarr(12,12,12)
lies=intarr(12,12,12)
options=intarr(12,12,12)
al=fltarr(12,12,12)
be=fltarr(12,12,12)
error=fltarr(12,12,12)

for k=0,1,1 do begin    ;划分R块   0-11
   for l=0,1,1 do begin   ;划分R块
     for h=0,1,1 do begin

      Rklh=long(xll[k*4:k*4+3,l*4:l*4+3,h*4:h*4+3])   ;划分R块
   ;   print,rklh
      rsum=total(Rklh)
      rsum2=total(Rklh*Rklh)


;下面对固定的一块Rklh，划分所有的Dijt与开匹配
        best_err_i=50000.0    ;用于保存最优误差值所在的Dij行
   		best_ops=0       ;用于保存最好0~23操作
   		best_Di=0        ; 从0行开始

     for i=0,40,1 do begin          ;0-44
        best_err_j=50000.0   ; 用于保存最优误差值所在的Dij列
        best_Dj=0

        for j=0,40,1 do begin  ;0-32
        best_err_t=50000.0
        best_Dt=0
           for t=0,40,1 do begin ;划分D块 ;0-44
             Dijt=xll[i:i+7,j:j+7,t:t+7]
            ;计算缩小的D块
            for m=0,3,1 do begin
              for n=0,3,1 do begin
                for q=0,3,1 do begin
                 DSijt(m,n,q)=mean(Dijt[2*m:2*m+1,2*n:2*n+1,2*q:2*q+1])
                endfor
              endfor
            endfor

;已完成D块的缩小，下面计算Rsum,rsum2,dsum,dsum2,rdsum

 dsum=total(DSijt)
 ;print,"dsum",dsum
 dsum2=total(long(DSijt)*long(DSijt))
 ;print,"dsum2",dsum2

;以下分别计算22种变换的rdsum及最佳参数
  	best_err_p=5000000.0    ;定义一个最优误差值
 	for p=0,21,1 do begin    ;24种变换
   		rdsum=0L
   		DSijt1=fltarr(4,4,4)
   		rd6sum2,Rklh,DSijt,DSijt1,rdsum,p
   		;print,"rdsum",rdsum
   	  ; 已返回第p种变换的rdsum
      ;计算alpha的分母
       ;  det=512*dsum2-dsum*dsum
       ;  if(det=0) then begin
        ;          alpha=0
        ; endif else begin

                    alpha=(64*rdsum-rsum*dsum)/(64*dsum2-dsum*dsum)
                  ;  print,"alhpa",alpha
           ; print,alpha
        ; endelse
             if (alpha lt 0) then alpha=0.05       ;    0<alpha<1
             if (alpha gt 1) then alpha=0.95    ;    0<alpha<1
             beta=(rsum-alpha*dsum)/64
            ; print,"beta",beta
             ;计算误差
            ; b=Rklh-alpha*DSijt-beta
             ;err=total(b*b)/512
             ;这里应该用变化后的DSijt1计算err
            err=total((Rklh-alpha*DSijt1-beta)*(Rklh-alpha*DSijt1-beta))/64
        ;    print,alpha
         ;   print,beta
          ;  print,err
;             rms= sqrt((rsum2 + alpha*(alpha*dsum2 - 2.0*rdsum + 2.0*beta*dsum) +
 ;                beta*(beta*w2 - 2.0*rsum))/w2);

           ; print,"err",err
             if(err lt best_err_p) then begin   ;在一块DSijt,  0-21种变换中寻找最优误差
                 best_err_p=err
                 best_alpha_p=alpha
                 best_beta_p=beta
                 best_ops_p=p
             endif
         ;    if (best_err_p lt 50000) then break     ;退出循环的条件

        endfor  ;for p

        if(best_err_p lt best_err_t) then begin  ;在一块DSijt， 0-21种变换后层中寻找最优误差
                 best_err_t=best_err_p
                 best_alpha_t=best_alpha_p
                 best_beta_t=best_beta_p
                 best_ops_t=best_ops_p
                 best_Dt=t        ;第几层
           endif
        ;   if (best_err_t lt 1300) then break    ;退出循环的条件

         endfor   ;for  t
        if(best_err_t lt best_err_j) then begin   ;在一块DSij,  0-21种变换后列中寻找最优误差
                 best_err_j=best_err_t
                 best_alpha_j=best_alpha_t
                 best_beta_j=best_beta_t
                 best_ops_j=best_ops_t
                 best_Dt_j=best_Dt    ;第j列中第t层
                 best_DJ=j        ;第j列
             endif
      ;    if (best_err_j lt 300) then break   ;退出循环的条件
         endfor   ;for j
             if(best_err_j lt best_err_i) then begin   ;在一块DSij,  0-21种变换后行中寻找最优误差
                 best_err_i=best_err_j
                 best_alpha_i=best_alpha_j
                 best_beta_i=best_beta_j
                 best_ops_i=best_ops_j
                 best_Dt_j_i=best_Dt_j   ;第i行j列中第t层
                 best_Dj_i=best_Dj       ;第i行中第j列
                 best_Di=i               ;第i行
            endif

      ;  if (best_err_i lt 300) then break

      endfor    ;for i
        ; print,"best_Di,best_Dj_i,best_Dt_j_i",best_Di,best_Dj_i,best_Dt_j_i
         ;print,"ops",best_ops_i
        ; print,"best_beta_i",best_beta_i
        ; print,"best_alpha_i",best_alpha_i
        ; print,"best_err_i",best_err_i

    cens(k,l,h)=best_Dt_j_i
    lies(k,l,h)=best_Dj_i
    hans(k,l,h)=best_Di
    options(k,l,h)=best_ops_i
    al(k,l,h)=best_alpha_i
    be(k,l,h)=best_beta_i
    error(k,l,h)=best_err_i
    endfor   ;for h

    endfor   ;for l
endfor    ;for k

openw,lun,'c:\cens.txt',/get_lun
printf,lun,cens,format='(11(i2,","),i2)'
openw,lun,'c:\lies.txt',/get_lun
printf,lun,lies,format='(11(i2,","),i2)'
openw,lun,'c:\hans.txt',/get_lun
printf,lun,hans,format='(11(i2,","),i2)'
openw,lun,'c:\options.txt',/get_lun
printf,lun,options,format='(11(i2,","),i2)'
openw,lun,'c:\alpha.txt',/get_lun
printf,lun,al,format='(11(f14,","),f14)'
openw,lun,'c:\bate.txt',/get_lun
printf,lun,be,format='(11(f14,","),f14)'
err=float(err)
openw,lun,'c:\err.txt',/get_lun
printf,lun,error,format='(11(f14,","),f14)'
print,"cens",cens
print,"lies",lies
print,"hans",hans
print,"options",options
print,"al",al
print,"be",be
print,"error",error
close,/all
t1=systime(1)
print,(t1-t0)/60,"min"
end