IDL实现非监督分类——k-mean

       原理：随机初始化五个聚类中心，几个都行，不要纠结怎么选，随便就行，但要拉开每个点的差距。然后遍历每个图像像素点，计算每个像素点到五个聚类中心的距离，离哪个近就将它归类到哪个类别里，遍历结束后，计算每类的平均位置，以平均位置为新的聚类中心再次聚类，聚类次数足够时，聚类中心将不再变化，图像完成分类。

       当然最好设置迭代次数，以防止聚类次数过多而导致运算时间过长。

pro kmean

filepath0='E:\yaoganhomework\ch03+IDL图像统计程序\ch03 IDL图像统计程序\jz512-01.TIF'
filepath1='E:\yaoganhomework\ch03+IDL图像统计程序\ch03 IDL图像统计程序\jz512-02.TIF'
filepath2='E:\yaoganhomework\ch03+IDL图像统计程序\ch03 IDL图像统计程序\jz512-03.TIF'
filepath3='E:\yaoganhomework\ch03+IDL图像统计程序\ch03 IDL图像统计程序\jz512-04.TIF'
filepath4='E:\yaoganhomework\ch03+IDL图像统计程序\ch03 IDL图像统计程序\jz512-05.TIF'
filepath5='E:\yaoganhomework\ch03+IDL图像统计程序\ch03 IDL图像统计程序\jz512-06.TIF'

jzimage=nonliner_strech(read_tiff(filepath0))
dims=size(jzimage, /DIMENSIONS)

jzimage1=intarr(dims[0],dims[1],6)
jzimage1[*,*,0]=nonliner_strech(read_tiff(filepath0))
jzimage1[*,*,1]=nonliner_strech(read_tiff(filepath1))
jzimage1[*,*,2]=nonliner_strech(read_tiff(filepath2))
jzimage1[*,*,3]=nonliner_strech(read_tiff(filepath3))
jzimage1[*,*,4]=nonliner_strech(read_tiff(filepath4))
jzimage1[*,*,5]=nonliner_strech(read_tiff(filepath5))

Var=intarr(6,1)
men=intarr(6,1)
halfvar=intarr(6,1)
twomen=intarr(6,1)
threevar=intarr(6,1)
;监督分类的各地物均值
;var=[83.52,37.14,34.52,82.32,54.22,22.16]
;men=[83.9,38.88,35.72,84.06,79.68,32.88]
;halfvar=[117.24,57.12,66.28,65,94.88,62.38]
;twomen=[110.54,54.84,50.52,28.84,24.12,13.92]
;threevar=[107.42,55.8,69.54,84.08,122.7,70.58]
for i=0,5 do begin
  ;方差
  var[i]=stddev(jzimage1[*,*,i])
  ;平均值
  men[i]=mean(jzimage1[*,*,i])
endfor
halfvar=0.5*var
twomen=2*men
threevar=3*var
print,var
print,men
print,halfvar
print,twomen
print,threevar

;定义迭代次数
dd=21
diedai=0

while diedai lt dd do begin
  fist=intarr(5,1)
  varcount=0
  mencount=0
  halfvarcount=0
  twomencount=0
  threevarcount=0
  v_ar=intarr(6,1)
  m_en=intarr(6,1)
  h_alfvar=intarr(6,1)
  t_women=intarr(6,1)
  t_hreevar=intarr(6,1)
  for i=0,dims[0]-1 do begin
    for j=0,dims[1]-1 do begin
      jz1=intarr(1,1,6)
      jz1[0,0,0]=jzimage1[i,j,0]
      jz1[0,0,1]=jzimage1[i,j,1]
      jz1[0,0,2]=jzimage1[i,j,2]
      jz1[0,0,3]=jzimage1[i,j,3]
      jz1[0,0,4]=jzimage1[i,j,4]
      jz1[0,0,5]=jzimage1[i,j,5]
      fist[0]=pf(var,jz1)
      fist[1]=pf(men,jz1)
      fist[2]=pf(halfvar,jz1)
      fist[3]=pf(twomen,jz1)
      fist[4]=pf(threevar,jz1)
      a=min(fist,index)
      case index of
        0:begin          
          v_ar[0]+=jzimage1[i,j,0]
          v_ar[1]+=jzimage1[i,j,1]
          v_ar[2]+=jzimage1[i,j,2]
          v_ar[3]+=jzimage1[i,j,3]
          v_ar[4]+=jzimage1[i,j,4]
          v_ar[5]+=jzimage1[i,j,5]
          varcount+=1
        end
        1:begin
          m_en[0]+=jzimage1[i,j,0]
          m_en[1]+=jzimage1[i,j,1]
          m_en[2]+=jzimage1[i,j,2]
          m_en[3]+=jzimage1[i,j,3]
          m_en[4]+=jzimage1[i,j,4]
          m_en[5]+=jzimage1[i,j,5]
          mencount+=1
        end
        2:begin
          h_alfvar[0]+=jzimage1[i,j,0]
          h_alfvar[1]+=jzimage1[i,j,1]
          h_alfvar[2]+=jzimage1[i,j,2]
          h_alfvar[3]+=jzimage1[i,j,3]
          h_alfvar[4]+=jzimage1[i,j,4]
          h_alfvar[5]+=jzimage1[i,j,5]
          halfvarcount+=1
        end
        3:begin
          t_women[0]+=jzimage1[i,j,0]
          t_women[1]+=jzimage1[i,j,1]
          t_women[2]+=jzimage1[i,j,2]
          t_women[3]+=jzimage1[i,j,3]
          t_women[4]+=jzimage1[i,j,4]
          t_women[5]+=jzimage1[i,j,5]
          twomencount+=1
        end
        4:begin
          t_hreevar[0]+=jzimage1[i,j,0]
          t_hreevar[1]+=jzimage1[i,j,1]
          t_hreevar[2]+=jzimage1[i,j,2]
          t_hreevar[3]+=jzimage1[i,j,3]
          t_hreevar[4]+=jzimage1[i,j,4]
          t_hreevar[5]+=jzimage1[i,j,5]
          threevarcount+=1
        end
      endcase
    endfor
  endfor
  diedai+=1
  ;重新计算聚类中心
  var=v_ar/varcount
  men=m_en/mencount
  halfvar=h_alfvar/halfvarcount
  twomen=t_women/twomencount
  threevar=t_hreevar/threevarcount
endwhile
for i=0,dims[0]-1 do begin
  for j=0,dims[1]-1 do begin
    jz1=intarr(1,1,6)
    jz1[0,0,0]=jzimage1[i,j,0]
    jz1[0,0,1]=jzimage1[i,j,1]
    jz1[0,0,2]=jzimage1[i,j,2]
    jz1[0,0,3]=jzimage1[i,j,3]
    jz1[0,0,4]=jzimage1[i,j,4]
    jz1[0,0,5]=jzimage1[i,j,5]
    fist[0]=pf(var,jz1)
    fist[1]=pf(men,jz1)
    fist[2]=pf(halfvar,jz1)
    fist[3]=pf(twomen,jz1)
    fist[4]=pf(threevar,jz1)
    a=min(fist,index)
    case index of
      0:begin
        jzimage1[i,j,0:2]=[255,0,0]
      end
      1:begin
        jzimage1[i,j,0:2]=[0,255,0]
      end
      2:begin
        jzimage1[i,j,0:2]=[0,0,255]
      end
      3:begin
        jzimage1[i,j,0:2]=[0,0,0]
      end
      4:begin
        jzimage1[i,j,0:2]=[255,255,255]
      end
    endcase
  endfor
endfor
b=image(jzimage1[*,*,0:2])
end

Function nonliner_strech, img_data

  dims = size(img_data, /DIMENSIONS)
  hist_array = HISTOGRAM(img_data,Min=0,Max=255)

  hist_array = hist_array*1.0/(dims[0]*dims[1])

  cumul_hist = fltarr(256)
  cumul_hist[0] = hist_array[0]

  For i=1,255 do begin

    cumul_hist[i] = hist_array[i]+cumul_hist [i-1]

  Endfor

  givemin=where(cumul_hist lt 0.02)
  givemax=where(cumul_hist gt 0.98)

  per_min = max(givemin)
  per_max = min(givemax)

  img_min = min(img_data)
  img_max = max(img_data)

  img_strech=round((img_data-per_min)*(245-10)/(per_max-per_min)+10.0)* $
    ((img_data ge per_min) and (img_data le per_max)) + $
    round((img_data-img_min)*(9-0)/(per_min-img_min)+0.0)* (img_data lt per_min) +$
    round((img_data-per_max)*(255-246)/(img_max-per_max)+246.0)* (img_data gt per_max)

  return,img_strech

  ;end img_var
End

function pf,point,jz1
p=0.0
p=((jz1[0,0,0]-point[0])*(jz1[0,0,0]-point[0])+(jz1[0,0,1]-point[1])*(jz1[0,0,1]-point[1])+$
  (jz1[0,0,2]-point[2])*(jz1[0,0,2]-point[2])+(jz1[0,0,3]-point[3])*(jz1[0,0,3]-point[3])+$
    (jz1[0,0,4]-point[4])*(jz1[0,0,4]-point[4])+(jz1[0,0,5]-point[5])*(jz1[0,0,5]-point[5]))
return,p
end

迭代五次

迭代十八次标题

迭代二十一次

 

 从效果来看，迭代次数与分类效果并没有直接关系，具体需要自行调试。