基于误差分析的组合分类器

asd

基于误差分析的组合分类器

陈学泓，陈晋，杨伟，朱锴.　基于误差分析的组合分类器研究[J]. 遥感学报, 2008, 12(5): 681-691

　提出了一种基于误差分析的组合分类器,通过结合两种监督分类方法,提出的算法分别估计了两
种监督分类方法在计算过程中的误差,给出了规则输出的置信区间,再根据置信区间的大小对两种分类方法
的输出结果进行加权平均,从而得到更精确的规则输出。利用该方法对遥感图像进行分类实验,在不同训练
样本分布与不同训练样本数量的情况下,比较新的组合分类器与单一分类器的精度。结果表明新的组合分
类器能够取得比单一的分类器更高的分类精度。结果还显示出,两个分类器的独立性越强,组合分类器的效
果越好。另外一个实验比较了新的组合分类器与和式规则组合分类器的分类精度,结果仍显示出了新方法
的优越性。



欢迎同行们测试，还望不吝指教。



程序代码(需在ENVI 4.3及以上版本中执行)，代码地址如下:

http://www.damipan.com/file/1HOzjvr.html

下载后，按照以下分类流程分类：

１.编译代码

２.打开影像，选取ＲＯＩ（类别数至少为4类）

３.运行代码

４.分类结果存如内存


或者复制以下代码

;------------------------------------------------------------------------------
;          Combined classifier program (Chen Xuehong)
;      writen by Chen Xuehong, beijing normal university, 2008-08
;
;The method combines the rule outputs of Support Vector Machine (In ENVI 4.3)
;and Maximumlikelihood classifier depending on the estimated error weights. The
;detail information can be referenced in:
;
;Chen Xuehong, Chen Jin, Yang Wei, Zhu Kai. Study on Combined Classifier Based
;on Error Analysis. Journal of Remote Sensing, 2008,12(5)
;------------------------------------------------------------------------------




PRO newclass, event

; Confidence level is set as 0.317
zhixindu=0.317


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;Read file;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Selectin Image file and Regions of Interest
envi_select, title='Input Filename', fid=fid, pos=pos
if (fid eq -1) then return
roi_ids = envi_get_roi_ids(fid=fid, roi_colors=roi_colors,roi_names=roi_names)
if (roi_ids[0] eq -1) then begin
    print, 'No regions associated with the selected file'
    return
endif
; Compound widget for ROI selection
base = widget_auto_base(title='ROI Selection')
wm   = widget_multi(base, list=roi_names, uvalue='list', /auto)
result = auto_wid_mng(base)
if (result.accept eq 0) then return
ptr = where(result.list eq 1, count)

; class name is set as the ROI_name
class_names = ['Unclassified', roi_names]
   num_classes = n_elements(roi_ids)

; At least four classes should be selected because of the principle of the method
if num_classes le 3 then begin
print,'Error: At least four classes should be selected'
return
endif

; Displaying color of the class is set as the color of ROI
lookup = bytarr(3,num_classes+1)
lookup[0,1] = roi_colors

; Read the Image file
envi_file_query, fid, fname=fname, ns=ns, nl=nl, nb=nb
dims = [-1,0,ns - 1 ,0,nl - 1]
ImgData=fltarr(ns,nl,nb)
For i = 0,nb-1 Do Begin
      Dt = Envi_Get_Data(Fid = Fid,dims = dims,pos=i)
      ImgData[*,*,i] = Dt[*,*]
EndFor
    map_info = envi_get_map_info(fid=fid)
pos=indgen(nb)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;Finish Reading file;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;



;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;Starting computing the combined rule outpust;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

; Output the rules of two-class SVM
ids=long(2)
ruleid=intarr(count,count)
svmr=fltarr(count,count,ns,nl) ;Maxtrix used to store the rule outputs of two-class SVM

for i=0, count-2 do begin
for j=i+1,count-1 do begin
   ids=[roi_ids[ptr],roi_ids[ptr[j]]]

; Two-class SVM in ENVI 4.3 is imported
   envi_doit, 'envi_svm_doit', $
      fid=fid, pos=pos, dims=dims, $
      /in_memory,R_FID=r_fid, $
      roi_ids=ids, RULE_FID=rule_fid, /RULE_IN_MEMORY,thresh=thresh, $
    penalty=penalty, kernel_type= kernel_type, $
      kernel_degree=kernel_degree, kernel_bias=kernel_bias
   ruleid[i,j]=rule_fid
   envi_file_query, ruleid[i,j], ns=ns, nl=nl, nb=nb
   dims = [-1,0,ns - 1 ,0,nl - 1]

; Read the rule output of Two-class SVM
   ImgData2=fltarr(ns,nl,nb)
      Dt1 = Envi_Get_Data(Fid = ruleid[i,j],dims = dims,pos=0)
        Dt2 = Envi_Get_Data(Fid = ruleid[i,j],dims = dims,pos=1)
   if roi_names le roi_names[j] then begin
        svmr[i,j,*,*] = Dt1[*,*]
        svmr[j,i,*,*] = Dt2[*,*]
    endif else begin
        svmr[i,j,*,*] = Dt2[*,*]
         svmr[j,i,*,*] = Dt1[*,*]
   endelse
   ENVI_FILE_MNG, ID=rule_fid, /REMOVE
   ENVI_FILE_MNG, ID=r_fid, /REMOVE
endfor
endfor

; Add the ENVI status reporting bar1
rstr = ["MULTI-CLASS SVM CLASSFIER (1 PROCESS LEFT)"]
envi_report_init, rstr, title=&quotROCESS 1 ", base=base
envi_report_inc, base, ns

; Compute the rule outputs and the estimated error of multi-class SVM
svmerr=fltarr(ns,nl,count)
psvm=fltarr(ns,nl,count)
psvmln=fltarr(ns,nl,count)
svmerror=fltarr(ns,nl,count,count)
countnum=count*(count-1)/2
f=F_CVF(zhixindu,count,countnum-count-1)

;Solve the overdetermined equations to get the rule outputs of multi-class SVM
for i=0,ns-1 do begin
   envi_report_stat,base, i, ns
for j=0,nl-1 do begin
   A=fltarr(count+1,count+1)

   for ai=0,count-1 do begin
   for aj=0,count-1 do begin
    if (ai eq aj) then A[ai,aj]= total(svmr[*,aj,i,j]^2)
    if (ai ne aj) then A[ai,aj]=-svmr[ai,aj,i,j]*svmr[aj,ai,i,j]
   endfor
   endfor

   A[0:count-1,count]=1.0
   A[count,0:count-1]=1.0
   l=fltarr(count+1)
   l[count]=1.0
   t=invert(A)
   ptemple=double(t) ## double(l)
   psvm[i,j,*]=ptemple[0:count-1]
   psvmln[i,j,*]=alog(psvm[i,j,*])   ;rule outputs of multi-class SVM?

;Compute the error of multi-class SVM
   z=fltarr(count,countnum)
   numad=1
   for zj=0,countnum-1 do begin
    zm = fix(zj+numad)/fix(count)
    zn = fix(zj+numad) mod fix(count)
    if (zn eq 0) then begin
     numad=numad+zm+1
     zm = fix(zj+numad)/fix(count)
     zn = fix(zj+numad) mod fix(count)
    endif
    z[zn,zj]=svmr[zm,zn,i,j]
    z[zm,zj]=-svmr[zn,zm,i,j]
   endfor

   Z=double(Z)
   s=total((Z ## ptemple[0:count-1])^2)/(countnum-count-1)
   svmerror[i,j,*,*]=s*invert( transpose (Z) ##Z)*f

   for k=0,count-1 do begin
    svmerr[i,j,k]=svmerror[i,j,k,k] ;Read the diagonal element of the error matrix
   endfor
    svmerr[i,j,*]=svmerr[i,j,*]/(psvm[i,j,*]^2) ;error of multi-class SVM
endfor
endfor


; Classified result of multi-class SVM
svmresult=lonarr(ns,nl)
for i=0,ns-1 do begin
for j=0,nl-1 do begin
   svmresult[i,j]=where (psvm[i,j,*] eq max(psvm[i,j,*]))
endfor
endfor


; Close the status bar1
envi_report_init, base=base, /finish


; Compute MLC classifier

; Read the Image file
envi_file_query, fid, ns=ns, nl=nl, nb=nb
dims = [-1,0,ns - 1 ,0,nl - 1]
ImgData=fltarr(ns,nl,nb)
For i = 0,nb-1 Do Begin
      Dt = Envi_Get_Data(Fid = Fid,dims = dims,pos=i)
      ImgData[*,*,i] = Dt[*,*]
   EndFor
    map_info = envi_get_map_info(fid=fid)

pos=indgen(nb)
cov=dblarr(nb,nb,count)
   meanvalue=dblarr(nb,count)
meanvalueerr=dblarr(nb,nb,count)
    pos=indgen(nb)

; Compute the statistics parameters of the ROIs(training sample) and the MCL error

; Read ROI data
for i=0L, count-1 do begin
    For k = 0,nb-1 Do Begin
        data = envi_get_roi_data(roi_ids[ptr], fid=fid, $
          pos=pos)
   endfor

    cov[*,*,i]=correlate(data,/covariance,/double) ; covariance matrix of roi

      for j=0,nb-1 do begin
      meanvalue[j,i]=mean(data[j,*])
    endfor

      n=size(data[j-1,*],/n_element)
      meanvalueerr[*,*,i]=cov[*,*,i]*nb*(n-1)/(n*(n-nb))*F_CVF(zhixindu,nb,n-nb)
endfor

meanvalue=double(meanvalue)
meanvalueerr=double(meanvalueerr)


probability=dblarr(ns,nl,count)
probabilitye=dblarr(ns,nl,count)
probabilityerr=dblarr(ns,nl,count)

; Add the status bar 2
rstr = ["MLC CLASSFIER (0 PROCESS LEFT)"]
envi_report_init, rstr, title="PROCESS 2", base=base
envi_report_inc, base, ns

; Compute the rule outputs and error of MLC
for i=0,ns-1 do begin
envi_report_stat,base, i, ns
for j=0,nl-1 do begin

   for k=0,count-1 do begin
    matrixtemp=reform(cov[*,*,k])
    arraytemp=dblarr(nb)
    for z=0,nb-1 do begin
     arraytemp[z]=Imgdata[i,j,z]-meanvalue[z,k]
    endfor
    arraytemp=double(arraytemp)
    probability[i,j,k]=alog(double(1.0/count))- 0.5*alog(determ(matrixtemp))-0.5*((arraytemp) ##invert(matrixtemp) ##transpose(arraytemp))-0.5*nb*alog(2*3.14)
    matrixtemp=invert(matrixtemp)
    errtemp=reform(meanvalueerr[*,*,k])
    probabilitye[i,j,k]=exp(probability[i,j,k])
    probabilityerr[i,j,k]=transpose(matrixtemp ## transpose(arraytemp)) ## errtemp ## (matrixtemp ## transpose(arraytemp))
   endfor

   sumtemp=total(probabilitye[i,j,*])
   probabilitye[i,j,*]=probabilitye[i,j,*]/sumtemp
   probability[i,j,*]=probability[i,j,*]-alog(sumtemp)    ;rule outputs of MLC
   probabilityerr[i,j,*]=float(probabilityerr[i,j,*])/(1+probabilitye[i,j,*])^2 ;Error of MLC
endfor
endfor

; Close the status bar 2
envi_report_init, base=base, /finish

probability=float(probability)
maxerr=fltarr(ns,nl,count)


; Compute the combined rule output

maxerr=probabilityerr*2
   ; As only mean value error is estimated and the covariance error is ignored in Chen X H (2008)
   ; An experienced paramter 2 is set to compensate the covariance error, which can make a better result


pfinal=fltarr(ns,nl,count)
pfinal= (maxerr * psvmln + svmerr * probability)/( maxerr+svmerr)

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;Finish computing the combined rule outputs;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;Output the classified result;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

classresult=intarr(ns,nl)
for i=0,ns-1 do begin
for j=0,nl-1 do begin
classresult[i,j]=where (pfinal[i,j,*] eq max(pfinal[i,j,*]))
endfor
endfor
classresult=fix(classresult)
classresult=classresult+1
help,classresult

class_names = class_names
lookup = lookup
bnames = ['NEW Classification']
file_type = ENVI_FILE_TYPE('ENVI Classification')

; The classified result is stored in the memory
ENVI_ENTER_DATA, classresult, num_classes=num_classes, $
   class_names=class_names, lookup=lookup, $
   file_type=file_type, bnames=bnames,map_info=map_info

;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;All finished;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;


END

liuxl

学习一下