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

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

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

GLCM IDL代码

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./content/3811116.html

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