图像纹理特征的学习

一．灰度差分统计法

设（x，y）为图像中的一点，该点和它只有微小距离的点

的灰度差值为：



用两张图片举一个例子：一张图片为墙面的纹理图（wall.jpg），另一张图片为大理石面纹理图(stone.jpg)。分别计算这两张图片的平均值，对比度，熵。





wall.png stone.png

matlab代码如下：

clear all;

clc;

J=imread('wall.jpg'); %读入纹理图像，分别输入wall.jpg和stone.jpg两幅图进行对比

A=double(J);

[m,n]=size(A); %求A矩阵的大小，赋值给m n

B=A;

C=zeros(m,n); %新建全零矩阵C，以下求解归一化的灰度直方图

for i=1:m-1

for j=1:n-1

B(i,j)=A(i+1,j+1);

C(i,j)=(round(A(i,j)-B(i,j)));

%源程序为C(i,j)=abs(round(A(i,j)-B(i,j))); 但是我觉得由公式（1）的计算应该不加abs.

end

end

h=imhist(mat2gray(C))/(m*n);

mean=0;con=0;ent=0; % 均值mean、对比度con和熵ent初始值赋零

for i=1:256 %循环求解均值mean、对比度con和熵ent

mean=mean+(i*h(i))/256;

con=con+i*i*h(i);

if(h(i)>0)

ent=ent-h(i)*log2(h(i));

end

end

mean,con,ent



由表格中得到的纹理特征可以分析：可以看出，两张图片熵的差距很大，当纹理越粗糙，那么所对应的熵就越大，反之，则越小。

二．自相关函数法

matlab代码如下：

建立.m文件。

function [epsilon,eta,C]=zxcor(f,D,m,n)

%自相关函数zxcor(),f为读入的图像数据，D为偏移距离，【m，n】是图像的尺寸数据

for epsilon=1:D %循环求解图像f(i,j)与偏离值为D的像素之间的相关值

for eta=1:D

temp=0;

fp=0;

for x=1:m

for y=1:n

if(x+ epsilon -1)>m|(y+ eta -1)>n

f1=0;

else

f1=f(x,y)*f(x+ epsilon -1,y+ eta -1);

end

temp=f1+temp;

fp=f(x,y)*f(x,y)+fp;

end

end

f2(epsilon, eta)=temp;

f3(epsilon, eta)=fp;

C(epsilon, eta)= f2(epsilon, eta)/ f3(epsilon, eta); %相关值C

end

end

epsilon =0:(D-1); %x方向的取值范围

eta =0:(D-1); %y方向的取值范围

下面继续用wall.jpg和stone.jpg两张图片进行相关函数法的分析。

matlab代码如下：

clear;

clc;

close all;

f11=imread('wall.jpg'); %读入砖墙面图像，图像数据赋值给f

f1=rgb2gray(f11); %彩色图像转换成灰度图像

f1=double(f1); %图像数据变为double类型

[m,n]=size(f1); %图像大小赋值为[m,n]

D=20; %偏移量为20

[epsilon1,eta1,C1]=zxcor(f1,D,m,n); %调用自相关函数

f22=imread('stone.jpg'); %读入大理石图像，图像数据赋值给f

f2=rgb2gray(f22);

f2=double(f2);

[m,n]=size(f2);

[epsilon2,eta2,C2]=zxcor(f2,20,m,n); %调用自相关函数

set(0,'defaultFigurePosition',[100,100,1000,500]); %修改图形图像位置的默认设置

set(0,'defaultFigureColor',[1 1 1]);

figure;subplot(121);imshow(f11);

subplot(122);imshow(f22);

figure;subplot(121);mesh(epsilon1,eta1,C1); %显示自相关函数与x，y的三维图像

xlabel(' epsilon ');ylabel(' eta '); %标示坐标轴变量

subplot(122);mesh(epsilon2,eta2,C2);

xlabel(' epsilon ');ylabel(' eta ');

二．灰度共生矩阵

Matlab的代码如下：

clear all;

clc;

Gray = imread('stone.jpg');% stone.jpg wall.jpg

[M,N,O] = size(Gray);

%1.将各颜色分量转化为灰度

%--------------------------------------------------------------------------

% Gray = double(0.3*Image(:,:,1)+0.59*Image(:,:,2)+0.11*Image(:,:,3));

%--------------------------------------------------------------------------

%2.为了减少计算量，对原始图像灰度级压缩，将Gray量化成16级

%--------------------------------------------------------------------------

for i = 1:M

for j = 1:N

for n = 1:256/16

if (n-1)*16<=Gray(i,j)&Gray(i,j)<=(n-1)*16+15

Gray(i,j) = n-1;

end

end

end

end

%--------------------------------------------------------------------------

%3.计算四个共生矩阵P,取距离为1，角度分别为0,45,90,135

%--------------------------------------------------------------------------

P = zeros(16,16,4);

for m = 1:16

for n = 1:16

for i = 1:M

for j = 1:N

if j<N&Gray(i,j)==m-1&Gray(i,j+1)==n-1

P(m,n,1) = P(m,n,1)+1;

P(n,m,1) = P(m,n,1);

end

if i>1&j<N&Gray(i,j)==m-1&Gray(i-1,j+1)==n-1

P(m,n,2) = P(m,n,2)+1;

P(n,m,2) = P(m,n,2);

end

if i<M&Gray(i,j)==m-1&Gray(i+1,j)==n-1

P(m,n,3) = P(m,n,3)+1;

P(n,m,3) = P(m,n,3);

end

if i<M&j<N&Gray(i,j)==m-1&Gray(i+1,j+1)==n-1

P(m,n,4) = P(m,n,4)+1;

P(n,m,4) = P(m,n,4);

end

end

end

if m==n

P(m,n,:) = P(m,n,:)*2;

end

end

end

%%---------------------------------------------------------

% 对共生矩阵归一化

%%---------------------------------------------------------

for n = 1:4

P(:,:,n) = P(:,:,n)/sum(sum(P(:,:,n)));

end

%--------------------------------------------------------------------------

%4.对共生矩阵计算能量、熵、惯性矩、相关4个纹理参数

%--------------------------------------------------------------------------

H = zeros(1,4);

I = H;

Ux = H; Uy = H;

deltaX= H; deltaY = H;

C =H;

for n = 1:4

E(n) = sum(sum(P(:,:,n).^2)); %%能量

for i = 1:16

for j = 1:16

if P(i,j,n)~=0

H(n) = -P(i,j,n)*log(P(i,j,n))+H(n); %%熵

end

I(n) = (i-j)^2*P(i,j,n)+I(n); %%惯性矩

Ux(n) = i*P(i,j,n)+Ux(n); %相关性中μx

Uy(n) = j*P(i,j,n)+Uy(n); %相关性中μy

end

end

end

for n = 1:4

for i = 1:16

for j = 1:16

deltaX(n) = (i-Ux(n))^2*P(i,j,n)+deltaX(n); %相关性中σx

deltaY(n) = (j-Uy(n))^2*P(i,j,n)+deltaY(n); %相关性中σy

C(n) = i*j*P(i,j,n)+C(n); %E(XY)

end

end

%%E(XY)-E(X)*E(Y) http://blog.csdn.net/lyflower/article/details/3104470
C(n) = (C(n)-Ux(n)*Uy(n))/sqrt(deltaX(n))/sqrt(deltaY(n)); %相关性 %%%% 0,45,90,135方向上的相关性依次为： 0.113814, 0.112050, 0.113968, 0.111966

end

%--------------------------------------------------------------------------

%求能量、熵、惯性矩、相关的均值和标准差作为最终8维纹理特征

%--------------------------------------------------------------------------

a1 = mean(E)

b1 = sqrt(cov(E))

a2 = mean(H)

b2 = sqrt(cov(H))

a3 = mean(I)

b3 = sqrt(cov(I))

a4 = mean(C)

b4 = sqrt(cov(C))

sprintf('0,45,90,135方向上的能量依次为： %f, %f, %f, %f',E(1),E(2),E(3),E(4)) % 输出数据;

sprintf('0,45,90,135方向上的熵依次为： %f, %f, %f, %f',H(1),H(2),H(3),H(4)) % 输出数据;

sprintf('0,45,90,135方向上的惯性矩依次为： %f, %f, %f, %f',I(1),I(2),I(3),I(4)) % 输出数据;

sprintf('0,45,90,135方向上的相关性依次为： %f, %f, %f, %f',C(1),C(2),C(3),C(4)) % 输出数据;

以stone.jpg为例运行的结果为：

0,45,90,135方向上的能量依次为： 0.572189, 0.548764, 0.562225, 0.545456

0,45,90,135方向上的熵依次为： 1.042052, 1.118325, 1.073047, 1.130627

0,45,90,135方向上的相关性依次为： 0.821187, 0.743442, 0.791676, 0.715337

0,45,90,135方向上的惯性矩依次为： 0.150355, 0.215821, 0.175188, 0.239465

以wall.jpg为例运行的结果为：

0,45,90,135方向上的能量依次为： 0.023909, 0.020111, 0.021309, 0.018415

0,45,90,135方向上的熵依次为： 4.044570, 4.247616, 4.188561, 4.326743

0,45,90,135方向上的相关性依次为： 0.845737, 0.766686, 0.797616, 0.720958

0,45,90,135方向上的惯性矩依次为： 2.137586, 3.234009, 2.805084, 3.867861

四．灰度-梯度共生矩阵

Matlab代码如下：

% 输入：IN-灰度图像

% 输出：OUT-纹理特征统计量矩阵

% 灰度梯度共生矩阵 H

%归一化灰度梯度矩阵 H_basic

%小梯度优势 T1

% 大梯度优势 T2

% 灰度分布的不均匀性 T3

% 梯度分布的不均匀性 T4

% 能量 T5

% 灰度平均 T6

% 梯度平均 T7

% 灰度均方差 T8

% 梯度均方差 T9

% 相关 T10

% 灰度熵 T11

% 梯度熵 T12

% 混合熵 T13

% 惯性 T14

% 逆差矩 T15

clc;

clear all;

IN=imread('D1_03.png');%D1_01.png

imshow(IN);

gray=256;

[R,C]=size(IN);

%采用平方求和计算梯度矩阵

GM=zeros(R-1,C-1);

for i=1:R-1

for j=1:C-1

n_GM=(IN(i,j+1)-IN(i,j))^2+(IN(i+1,j)-IN(i,j))^2;

GM(i,j)=sqrt(double(n_GM));

end

end

% figure,imshow(GM);

%找出最大值最小值

n_min=min(GM(:));

n_max=max(GM(:));

%把梯度图像灰度级离散化

%设置新的灰度级为new_gray

new_gray=32;

%新的梯度矩阵为new_GM

new_GM=zeros(R-1,C-1);

new_GM=uint8((GM-n_min)/(n_max-n_min)*(new_gray-1));

%计算灰度梯度共生矩阵

%梯度矩阵比灰度矩阵维数少1，忽略灰度矩阵最外围

H=zeros(gray,new_gray);

for i=1:R-1

for j=1:C-1

H(IN(i,j)+1,new_GM(i,j)+1)= H(IN(i,j)+1,new_GM(i,j)+1)+1;

end

end

%归一化灰度梯度矩阵 H_basic

total=i*j;

H_basic=H/total;

%小梯度优势 T1

TT=sum(H);

T1=0;

for j=1:new_gray

T1=T1+TT(1,j)/j^2;

end

T1=T1/total;

%计算大梯度优势 T2

T2=0;

for j=1:new_gray

T2=T2+TT(1,j)*(j-1);

end

T2=T2/total;

%计算灰度分布的不均匀性 T3

T3=0;

TT1=sum(H');

for j=1:gray

T3=T3+TT1(1,j)^2;

end

T3=T3/total;

%计算梯度分布的不均匀性 T4

T4=0;

for j=1:new_gray

T4=T4+TT(1,j)^2;

end

T4=T4/total;

%计算能量 T5

T5=0;

for i=1:gray

for j=1:new_gray

T5=T5+H_basic(i,j)^2;

end

end

%计算灰度平均 T6

TT2=sum((H_basic)');

T6=0;

for j=1:gray

T6=T6+(j-1)*TT2(1,j);

end

%计算梯度平均 T7

T7=0;

TT3=sum(H_basic);

for j=1:new_gray

T7=T7+(j-1)*TT3(1,j);

end

%计算灰度均方差 T8

T8=0;

for j=1:gray

T8=T8+(j-1-T6)^2*TT2(1,j);

end

T8=sqrt(T8);

%计算梯度均方差 T9

T9=0;

for j=1:new_gray

T9=T9+(j-1-T7)^2*TT3(1,j);

end

T9=sqrt(T9);

% 计算相关 T10

T10=0;

for i=1:gray

for j=1:new_gray

T10=T10+(i-1-T6)*(j-1-T7)*H_basic(i,j);

end

end

%计算灰度熵 T11

T11=0;

for j=1:gray

T11=T11+TT2(1,j)*log10(TT2(1,j)+eps);

end

T11=-T11;

%计算梯度熵 T12

T12=0;

for j=1:new_gray

T12=T12+TT3(1,j)*log10(TT3(1,j)+eps);

end

T12=-T12;

%计算混合熵 T13

T13=0;

for i=1:gray

for j=1:new_gray

T13=T13+H_basic(i,j)*log10(H_basic(i,j)+eps);

end

end

T13=-T13;

%计算惯性 T14

T14=0;

for i=1:gray

for j=1:new_gray

T14=T14+(i-j)^2*H_basic(i,j);

end

end

%计算逆差矩 T15

T15=0;

for i=1:gray

for j=1:new_gray

T15=T15+H_basic(i,j)/(1+(i-j)^2);

end

end

%x=1:50:750;

sprintf('小梯度优势T1为：%f',T1)

OUT(1,1)=T1;

sprintf('大梯度优势T2为：%f',T2)

OUT(1,2)=T2;

sprintf('灰度分布不均匀性T3为:%f',T3)

OUT(1,3)=T3;

sprintf('梯度分布不均匀性T4为：%f',T4)

OUT(1,4)=T4;

sprintf('能量T5为：%f',T5)

OUT(1,5)=T5;

sprintf('灰度平均T6为：%f',T6)

OUT(1,6)=T6;

sprintf('梯度平均T7为：%f',T7)

OUT(1,7)=T7;

sprintf('灰度方差T8为：%f',T8)

OUT(1,8)=T8;

sprintf('梯度方差T9为：%f',T9)

OUT(1,9)=T9;

sprintf('相关T10为：%f',T10)

OUT(1,10)=T10;

sprintf('灰度熵T11为：%f',T11)

OUT(1,11)=T11;

sprintf('梯度熵T12为：%f',T12)

OUT(1,12)=T12;

sprintf('混合熵T13为：%f',T13)

OUT(1,13)=T13;

sprintf('惯性T14为：%f',T14)

OUT(1,14)=T14;

sprintf('逆差距T15为：%f',T15)

OUT(1,15)=T15;

对Brodatz纹理图像库1.1中的D1_01.png进行计算得到T1到T15的数据













从这几张图片中可以看出相同类型的纹理，所计算的纹理特征统计值（T1到T15）基本上差距不大，而纹理不同的所计算的值相差较大。

因此后期我打算将所有图片计算出来的T1到T15的值全部当做图片的特征向量，然后通过分类的方法来识别未知的纹理图片是属于哪一类的。