灰度共生矩阵GLCM及纹理特征影像生成

灰度共生矩阵GLCM及纹理特征影像生成

　　实现类似于滤波过程中的5*5窗体移动，形成子图像的过程，这里的方法边界的象元，滑动窗口元素补0：

Pro Texture
;针对灰度影像
file=Dialog_Pickfile(Filter='*.bmp',/Must_exist)
image=Read_Bmp(file)
sz=size(image)
m=sz[1]
n=sz[2]
;影像压缩成16个灰度级
dimage=uintarr(m,n)
for i=0,m-1 do begin
  for j=0,n-1 do begin
  dimage[i,j]=image[i,j] * 16 / 256
  endfor
endfor
;构造子影像，窗口大小5X5
subimage=uintarr(5,5)
e=0
for i=0,m-1 do begin
  for j=0,n-1 do begin
    ;填充窗口内的值
    for k=-2,2,1 do begin
      t=i+k
      if t Lt 0 then begin
          Continue
      endif
      if t ge m then begin
          Break
      Endif
      for l=-2,2,1 do begin
        if(j+l) lt 0 then begin
          continue;
        endif
        if(j+l) lt n then begin
          ;给子影像赋值
          subimage[k+2,l+2]=dimage[i+k,j+l];
        endif
      endfor
    endfor
;   test=uintarr(5,5);
;   for r=0,4 do begin
;   for s=0,4 do begin
;      test[r,s]=subimage[r,s]
;   Endfor
;   endfor
;应用子影像进行计算

  endfor
endfor
TV,image,True =1
end

　　下面是glcm计算的第一个版本，直接用IDL编程有两个问题需要注意，一个是逻辑运算符，IDL中的>和<符号不表示比较，而是表示取大和取消，逻辑运算采用字母表示，如ge，gt等；二是两个整数相除，得到的结果仍是整数，需要将被除数类型转换。

Function Glcm,subImage,d
sz=size(subImage)
m=sz[1]
n=sz[2]
glcmH=lonarr(16,16)
glcmEN=lonarr(16,16)
glcmV=lonarr(16,16)
glcmWN=lonarr(16,16)

;0度
for i= 0,m-1 do begin
  for j=0,n-d-1 do begin
     glcmH[subImage[i,j],subImage[i,j+d]]=glcmH[subImage[i,j],subImage[i,j+d]]+1;%是共生矩阵0度的计算式
     glcmH[subImage[i,j+d],subImage[i,j]]=glcmH[subImage[i,j+d],subImage[i,j]]+1
  endFor
endFor
;45度
for ii= 0,m-d-1 do begin
  for jj=0,n-d-1 do begin
     glcmEN[subImage[ii,jj],subImage[ii+d,jj+d]]=glcmEN[subImage[ii,jj],subImage[ii+d,jj+d]]+1;%是共生矩阵45度的计算式
     glcmEN[subImage[ii+d,jj+d],subImage[ii,jj]]=glcmEN[subImage[ii+d,jj+d],subImage[ii,jj]]+1
  endFor
endFor
;90度
for iii= 0,m-d-1 do begin
  for jjj=0,n-1 do begin
     glcmV[subImage[iii,jjj],subImage[iii+d,jjj]]=glcmV[subImage[iii,jjj],subImage[iii+d,jjj]]+1;%是共生矩阵90度的计算式
     glcmV[subImage[iii+d,jjj],subImage[iii,jjj]]=glcmV[subImage[iii+d,jjj],subImage[iii,jjj]]+1
  endFor
endFor
;135度
for iiii= d,m-1 do begin
  for jjjj=0,n-d-1 do begin
     glcmWN[subImage[iiii,jjjj],subImage[iiii-d,jjjj+d]]=glcmWN[subImage[iiii,jjjj],subImage[iiii-d,jjjj+d]]+1;%是共生矩阵135度的计算式
     glcmWN[subImage[iiii-d,jjjj+d],subImage[iiii,jjjj]]=glcmWN[subImage[iiii-d,jjjj+d],subImage[iiii,jjjj]]+1
  endFor
endFor
;归一化
;sumH=0L
;sumEN=0L
;sumV=0L
;sumWN=0L
;for k=0,15 do begin
;  for l=0,15 do begin
;  sumH=sumH+glcmH[k,l]
;  sumEN=sumEN+glcmEN[k,l]
;  sumV=sumV+glcmV[k,l]
;  sumWN=sumWN+glcmWN[k,l]
;  endfor
;endfor
glcmHU=dblarr(16,16)
glcmENU=dblarr(16,16)
glcmVU=dblarr(16,16)
glcmWNU=dblarr(16,16)
for kk=0,15 do begin
  for ll=0,15 do begin
  glcmHU[kk,ll]=double(glcmH[kk,ll])/(2*m*(n-d))
  glcmENU[kk,ll]=double(glcmEN[kk,ll])/(2*(m-d)*(n-d))
  glcmVU[kk,ll]=double(glcmV[kk,ll])/(2*(m-d)*n)
  glcmWNU[kk,ll]=double(glcmWN[kk,ll])/(2*(m-d)*(n-d))
  endfor
endfor
;求取特征值
e=dblarr(4)
H=dblarr(4)
In=dblarr(4)
  Ux=dblarr(4)
  Uy=dblarr(4)
  deltaX=dblarr(4)
  deltaY=dblarr(4)
C=dblarr(4)
for i1=0,15 do begin
  for j1=0,15 do begin
  e[0]=e[0]+glcmHU[i1,j1]*glcmHU[i1,j1]
  e[1]=e[1]+glcmENU[i1,j1]*glcmENU[i1,j1]
  e[2]=e[2]+glcmVU[i1,j1]*glcmVU[i1,j1]
  e[3]=e[3]+glcmWNU[i1,j1]*glcmWNU[i1,j1]
  
   if glcmHU[i1,j1] ne 0 then begin
      H[0] = -glcmHU[i1,j1]*Alog10(glcmHU[i1,j1])+H[0]; %%熵     
   endif
   if glcmENU[i1,j1] ne 0 then begin
      H[1] = -glcmENU[i1,j1]*Alog10(glcmENU[i1,j1])+H[1]; %%熵      
   endif
   if glcmVU[i1,j1] ne 0 then begin
      H[2] = -glcmVU[i1,j1]*Alog10(glcmVU[i1,j1])+H[2]; %%熵    
   endif
   if glcmWNU[i1,j1] ne 0 then begin
      H[3] = -glcmWNU[i1,j1]*Alog10(glcmWNU[i1,j1])+H[3]; %%熵     
   endif
   In[0] = (i1-j1)^2*glcmHU[i1,j1]+In[0];  %%惯性矩
   In[1] = (i1-j1)^2*glcmENU[i1,j1]+In[1];
   In[2] = (i1-j1)^2*glcmVU[i1,j1]+In[2];
   In[3] = (i1-j1)^2*glcmWNU[i1,j1]+In[3];
   
   Ux[0] = i1*glcmHU[i1,j1]+Ux[0]; %相关性中μx
   Uy[0] = j1*glcmHU[i1,j1]+Uy[0]; %相关性中μy
   Ux[1] = i1*glcmENU[i1,j1]+Ux[1]; %相关性中μx
   Uy[1] = j1*glcmENU[i1,j1]+Uy[1]; %相关性中μy
   Ux[2] = i1*glcmVU[i1,j1]+Ux[2]; %相关性中μx
   Uy[2] = j1*glcmVU[i1,j1]+Uy[2]; %相关性中μy
   Ux[3] = i1*glcmWNU[i1,j1]+Ux[3]; %相关性中μx
   Uy[3] = j1*glcmWNU[i1,j1]+Uy[3]; %相关性中μy
  endfor
endfor
for i2=0,15 do begin
  for j2=0,15 do begin
  deltaX[0] = (i2-Ux[0])^2*glcmHU[i2,j2]+deltaX[0]; %相关性中σx
  deltaY[0] = (j2-Uy[0])^2*glcmHU[i2,j2]+deltaY[0]; %相关性中σy
  C[0] = i2*j2*glcmHU[i2,j2]+C[0]; 
  
  deltaX[1] = (i2-Ux[1])^2*glcmENU[i2,j2]+deltaX[1]; %相关性中σx
  deltaY[1] = (j2-Uy[1])^2*glcmENU[i2,j2]+deltaY[1]; %相关性中σy
  C[1] = i2*j2*glcmENU[i2,j2]+C[1]; 
  
  deltaX[2] = (i2-Ux[2])^2*glcmVU[i2,j2]+deltaX[2]; %相关性中σx
  deltaY[2] = (j2-Uy[2])^2*glcmVU[i2,j2]+deltaY[2]; %相关性中σy
  C[2] = i2*j2*glcmVU[i2,j2]+C[2]; 
  
  deltaX[3] = (i2-Ux[3])^2*glcmWNU[i2,j2]+deltaX[3]; %相关性中σx
  deltaY[3] = (j2-Uy[3])^2*glcmWNU[i2,j2]+deltaY[3]; %相关性中σy
  C[3] = i2*j2*glcmWNU[i2,j2]+C[3]; 
  endfor
endfor
C[0] = (C[0]-Ux[0]*Uy[0])/deltaX[0]/deltaY[0]; %相关性
C[1] = (C[1]-Ux[1]*Uy[1])/deltaX[1]/deltaY[1];
C[2] = (C[2]-Ux[2]*Uy[2])/deltaX[2]/deltaY[2];
C[3] = (C[3]-Ux[3]*Uy[3])/deltaX[3]/deltaY[3];
;特征值求平均
feature=dblarr(4)
for i3=0,3 do begin
  feature[0]=feature[0]+e[i3]
  feature[1]=feature[1]+H[i3]
  feature[2]=feature[2]+In[i3]
  feature[3]=feature[3]+C[i3]
endfor
for i4=0,3 do begin
  feature[i4]=feature[i4]/4
endfor
return,feature;返回特征值
end

　　该函数传入灰度图像subImage和步长参数d，调用GLCM函数的测试代码：

pro test_glcm
file=Dialog_Pickfile(Filter='*.bmp',/Must_exist)
image=Read_Bmp(file)
sz=size(image)
m=sz[1]
n=sz[2]
;影像压缩成16个灰度级
dimage=uintarr(m,n)
textImage=uintarr(m,n)
for i=0,m-1 do begin
  for j=0,n-1 do begin
  dimage[i,j]=image[i,j] * 16 / 256
  endfor
endfor
feature=glcm(dimage,1)
print,feature
end

这里贴出代码，希望熟悉的朋友也帮忙验证一下，看看能否对其中的某些部分优化一下！

参考文献：

1.贾永红《数字图像处理》武汉大学出版社 P182~P184

2.http://www.cnblogs.com/skyseraph/archive/2011/08/27/2155776.html

3.http://blog.sina.com.cn/s/blog_4b700c4c0102e038.html