最大似然分类器代码

;最大似然分类器
;遥感11-3班 许忠雄 2014年6月1日（编）
pro MLC
Compile_opt idl2
;启动ENVI二次开发
ENVI,/restore_base_save_files
ENVI_BATCH_INIT
;文件绝对路径
tm_file ='E:\Classification_tm\LT51240362007139_JZ.img'
;读取文件
ENVI_OPEN_FILE, tm_file, r_fid=fid
;获取文件信息
ENVI_FILE_QUERY, fid, dims=dims, ns=ns, nl=nl, nb=nb
;创建影像数组，并将原始图像赋值进去
tm_data = intarr(ns,nl,nb)
tm_data[*,*,0] = envi_get_data(fid=fid, dims=dims, pos=0)
tm_data[*,*,1] = envi_get_data(fid=fid, dims=dims, pos=1)
tm_data[*,*,2] = envi_get_data(fid=fid, dims=dims, pos=2)
tm_data[*,*,3] = envi_get_data(fid=fid, dims=dims, pos=3)
tm_data[*,*,4] = envi_get_data(fid=fid, dims=dims, pos=4)
tm_data[*,*,5] = envi_get_data(fid=fid, dims=dims, pos=5)
;获取影像投影信息
map_info = envi_get_map_info(fid=fid)
;创建分类后的影像数组
class_data=intarr(ns,nl,3)
;训练样本，根据自己的需要选择一定数量的类别，本程序是4个类别
;分别为水体、裸地、植被和城市，每类在各波段中取8个点（为了提高精度可以将点数增多）
;以下数组列数等于波段数，行数等于选择点数
water=[[117,59,55,29,25,13]$;水体
,[115,59,54,29,24,14]$
,[103,52,46,26,24,14]$
,[117,58,55,29,25,14]$
,[111,54,49,28,23,13]$
,[107,53,48,27,22,12]$
,[111,55,51,29,24,12]$
,[107,53,50,28,22,12]]
bare_land=[[92,45,58,76,97,48]$;裸地
,[94,47,59,71,95,51]$
,[99,50,66,72,107,59]$
,[98,49,64,79,99,52]$
,[96,47,59,69,95,52]$
,[95,47,59,86,102,52]$
,[94,47,57,78,97,49]$
,[95,46,58,76,101,53]]
plant=[[85,39,37,92,64,27]$;植被
,[82,39,36,94,71,28]$
,[83,37,36,88,64,26]$
,[84,37,35,95,61,21]$
,[88,41,41,86,75,31]$
,[86,38,37,89,71,29]$
,[83,37,36,84,62,26]$
,[88,40,39,94,71,28]]
urban=[[108,48,53,56,81,47]$;城市
,[111,53,60,66,89,55]$
,[115,57,67,64,99,69]$
,[128,64,73,68,97,63]$
,[101,49,52,49,71,48]$
,[116,59,68,64,96,62]$
,[115,55,65,61,96,63]$
,[113,55,64,66,89,58]]

;………………（根据自己类的数量在这添加代码）

;创建进度条，算法比较复杂，运行时间长，创建进度条心里有个底--（纯属个人喜欢）
tlb=widget_base(xsize=330,ysize=100,xoffset=500,yoffset=250,title='耐心等待！')
;显示进度条
widget_control,tlb,/real
prsbar=idlitwdprogressbar(group_leader=tlb,title='分类中…',cancel=cancelin)
;求各类在各波段中的均值
;class_mean 为自己编写的函数，在这进行调用，函数是实现在最下面
mean_water=class_mean(water)
mean_bare_land=class_mean(bare_land)
mean_plant=class_mean(plant)
mean_urban=class_mean(urban)

;………………（根据自己类的数量在这添加代码）

;求各类的协方差矩阵，直接调用IDL内部函数imsl_covariances
cov_water=imsl_covariances(transpose(water))
cov_bare_land=imsl_covariances(tran

spose(bare_land))
cov_plant=imsl_covariances(transpose(plant))
cov_urb
an=imsl_covariances(transpose(urban))

;………………（根据自己类的数量在这添加代码）

;遍历图像所有行所有列
for i=0,ns-1 do begin
for j=0,nl-1 do begin
;根据最大似然分类器求个像元落在各类的“概率”（简化后的概率）
f_water=alog(determ(cov_water))+transpose(tm_data[i,j,*]-$
mean_water)##invert(cov_water)##transpose(reform(tm_data[i,j,*]-mean_water))

f_bare_land=alog(determ(cov_bare_land))+transpose(tm_data[i,j,*]-$
mean_bare_land)##invert(cov_bare_land)##transpose(reform(tm_data[i,j,*]-mean_bare_land))

f_plant=alog(determ(cov_plant))+transpose(tm_data[i,j,*]-$
mean_plant)##invert(cov_plant)##transpose(reform(tm_data[i,j,*]-mean_plant))

f_urban=alog(determ(cov_urban))+transpose(tm_data[i,j,*]-$
mean_urban)##invert(cov_urban)##transpose(reform(tm_data[i,j,*]-mean_urban))

;………………（根据自己类的数量在这添加代码）

pro_mtx=[f_water,f_bare_land,f_plant,f_urban];,……]
;由于是取简化后的“概率”式子，即仅取被减掉的部分，所以在这取最小值！
min_pro=min(pro_mtx)
;检索像元落在哪个类中并给其重新赋值R、G、B
if min_pro eq f_water then class_data[i,j,*]=[23,23,115]
if min_pro eq f_bare_land then class_data[i,j,*]=[239,183,102]
if min_pro eq f_plant then class_data[i,j,*]=[49,205,49]
if min_pro eq f_urban then class_data[i,j,*]=[255,0,0]

;………………（根据自己类的数量在这添加代码）

endfor
;进度条显示处理进度
idlitwdprogressbar_setvalue,prsbar,(i*1.0)*100/ns
endfor
;将处理后的图像写出来！
pos_class=indgen(3)
out_name ='E:\Classification_tm\tm_class_jz.img'
envi_enter_data, class_data, map_info = map_info, r_fid=out_fid
envi_output_to_external_format, fid=out_fid, out_name = out_name, dims=dims, pos=pos_class, /envi
widget_control,tlb,/destroy
Print, 'END all!!!'
end

;自己定义的求各类均值（求各列的均值）的函数，传递参数为类名
function class_mean,class_name
dims=size(class_name,/DIMENSIONS)
class_mean=fltarr(dims[0])
for i=0,dims[0]-1 do begin
class_mean[i]=mean(class_name[i,*])
endfor
return,class_mean
end