Tamura 纹理特征

我这篇文章主要是参考的Tamura纹理特征的matlab实现。本来没打算写这篇博客的，结果在写文章的时候各种找文献资料，都很难找到比较好的解释Tamura的文章。很多人的文章都是含糊其辞，要么就是排版稀烂，没法看。实在受不了自己写一个高大全的Tamura特征的博客，既方便自己也方便别人。

• 原始出处
• 原理解释
• 代码展示

原始出处

最原始的Tamura的论文《Textural Features Corresponding to Visual Perception》，IEEE上的这篇论文，可以直接下载。
Tamura纹理特征在一些中文论文中也有出现，也有解释的不错的，我推荐《Tamura纹理特征在水下目标分类中的应用》但是基本都避开了Tamura特征中线性相关的三个特征。
我在网上还找了一个讲纹理特征计算的博客，里面讲了Tamura纹理特征，还算是比较细致的，也不错。

原理解释

Tamura纹理特征包括六个指标：$粗糙度（Coarseness）、对比度（ Contrast）、方向度（ Directionality）、线性度（ Linelikeness）、规则度（Regularity）、粗略度（Roughness）$，一般论文里面只用前三个特征，说前面三个特征是线性无关的，后面三个特征和前面三个特征是线性相关的，因此只采用前三个特征。但是我发现其实线性度这个特征还是比较难解释的，单纯观察看到的代码而言。

$粗糙度（Coarseness）$

粗糙度是反映纹理中粒度的一个量，是最基本的纹理特征。当两种纹理特征模式知识基元尺寸不同时，具有较大基元尺寸的模式给人感觉更粗糙。粗糙度的计算可以分为以下几个步骤进行：
.
1. 计算图像中大小为$2^k×2^k$个像素的活动窗口中像素的平均强度值，公式：
$A_{k}(x,y)=\sum_{i=x-2^{k-1}}^{x+2^{k-1}-1}\sum_{j=y-2^{k-1}}^{y+2^{k-1}-1}g(i,j)/2^{2k}$

式中：$k=0,1,...,5;g(i,j)是位于(i,j)处像素的灰度值。$
.
2. 对每个像素，分别计算在水平和垂直方向上互不重叠的窗口之间的平均强度差，公式：
$E_{k,h}(x,y)=\begin{vmatrix} A_{k}(x+2^{k-1},y)-A_{k}(x-2^{k-1},y) \end{vmatrix}, E_{k,v}(x,y)=\begin{vmatrix} A_{k}(x,y+2^{k-1})-A_{k}(x,y-2^{k-1}) \end{vmatrix}$
其中对于每个像素，使得$E$值达到最大的$k$值用来设置最佳尺寸$S_{best}(x,y)=2^k$.
.
3. 粗糙度通过计算整幅图像中的$S_{best}$的平均值得到，公式：
$F_{crs}=\frac{1}{m×n}\sum_{m}^{i=1}\sum_{n}^{j=1}S_{best}(i,j)$

$对比度（ Contrast）$

对比度是通过对像素强度分布情况的统计得到的，其大小由四个因素决定：灰度动态范围、直方图上黑白部分两极分化程度、边缘锐度和重复模式的周期。一般情况下，对比度指前面两个因素。
公式：$F_{con}=\frac{\sigma }{\alpha _{1/4}^{4}}$
其中：$\alpha _{4}=\mu _{4}/\sigma ^{4},\mu _{4}是四次矩而\sigma ^{2}是方差$
对比度给出了整幅图像的全局度量

$方向度（ Directionality）$

方向度是给定纹理区域的全局特性，描述纹理如何沿着某些方向发散或者集中的。计算步骤：
1. 每个像素出的梯度向量
$\left | \Delta G \right |=(\left |\Delta H \right |+\left |\Delta V \right |)/2$
$\Theta =tan^{-1}(\Delta V /\Delta H)+\pi /2$
其中：$\Delta V 和\Delta H$分别是通过图像卷积下面两个$3×3$卷积核得到的水平和垂直方向上的变化量。

当所有的像素梯度向量都被计算出来后，可以使用直方图$H_D表达\theta值。$
2. 通过直方图对$\theta 值域范围进行离散化$
3. 统计每个柄中相应的$\left| \Delta G \right|$大于给定阈值的像素数量。 对于有明显方向性的图像会出现峰值，对于无明显方向性的图像表现的比较平坦。
4. 图像总体的方向性可以通过计算直方图中峰值的尖锐程度获得，公式：
$F_{dir}=\sum_{n_{p}}^{p}\sum_{\Phi \epsilon W_{p}}(\Phi -\Phi _{p})^{2}H_{D}(\Phi )$
上式中$p$代表直方图中的峰值，$n_p$为直方图中所有的峰值，对于某个峰值$p$，$W_p$代表该峰值所包含的所有的离散区域，而$\Phi_p$是波峰中心位置

$线性度（ Linelikeness）$

线性度几乎所有的中文文献上都没有提及，今天我这相当于独一份儿了。（接下来的内容都是我从原始论文上搬下来的，可能有些地方翻译的不到位，大家还望见谅）
似乎我们需要对线性度提出更多细节方面的描述。我们所指的线性度是由线段组成的纹理的一种情况。我们认为一个方向以及邻近的方向角度近似相等，我们认为这样的一组边缘点集是一条线段上的。
为了更加具体的描述，我们建立了方向共生矩阵，其中的元素$P_{Dd}(i,j)$是矩阵元素，它代表在图像中沿着边缘方向相距$d$的相邻两像素出现的频率，其中一个方向编码$i$，另一个的方向编码是$j$(本人认为这是代表不同的方向角度)。

这个在图中很明显的表示出来了。好几个纹理特征可以从中推导出来。我们定义如下的测度规则，共生矩阵中同一方向的用$+1$表示，垂直方向的用$-1$表示。
$F_{lin}=\sum_{i}^{n}\sum_{j}^{n}P_{Dd}(i,j)cos[(i-j)\frac{2\pi }{n}]/\sum_{i}^{n}\sum_{j}^{n}P_{Dd}(i,j)$
其中：$P_{Dd}表示在距离内的 n×n大小的局部方向共生矩阵$

$规则度（Regularity）$

公式：$F_{reg}=1-r(\sigma _{crs}+\sigma _{con}+\sigma _{dir}+\sigma _{lin})$
其中：$r$是归一化因子，每一个$\sigma_{xxx}都是代表F_{xxx}的标准差$

$粗略度（Roughness）$

公式：$F_{rgh}=F_{crs}+F_{con}$

代码展示

都是MATLAB写的函数代码

function feature=Tamura(imag)
             feature=zeros(1,5);   
             Fcrs=coarseness(imag,4);
             Fcon=contrast(imag);
             [Fdir,sita]=directionality(imag);
             Flin=linelikeness(imag,sita,4);
            %  Freg{k,1}=regularity(I{1,k},32);
             Frgh=Fcrs+Fcon;
             feature(1,:)=[ Fcrs, Fcon,Fdir,Flin,Frgh];
end

%% 第一个指标 Coarseness，粗糙度

function Fcrs = coarseness( graypic,kmax )%graphic为待处理的灰度图像，2^kmax为最大窗口  
[h,w]=size(graypic); %获取图片大小  
A=zeros(h,w,2^kmax); %平均灰度值矩阵A  
%计算有效可计算范围内每个点的2^k邻域内的平均灰度值  
for i=2^(kmax-1)+1:h-2^(kmax-1)  
    for j=2^(kmax-1)+1:w-2^(kmax-1)  
        for k=1:kmax  
            A(i,j,k)=mean2(graypic(i-2^(k-1):i+2^(k-1)-1,j-2^(k-1):j+2^(k-1)-1));  
        end  
    end  
end  
%对每个像素点计算在水平和垂直方向上不重叠窗口之间的Ak差  
for i=1+2^(kmax-1):h-2^(kmax-1)  
    for j=1+2^(kmax-1):w-2^(kmax-1)  
        for k=1:kmax  
            Eh(i,j,k)=abs(A(i+2^(k-1),j,k)-A(i-2^(k-1),j));  
            Ev(i,j,k)=abs(A(i,j+2^(k-1),k)-A(i,j-2^(k-1)));  
        end  
    end  
end  
%对每个像素点计算使E达到最大值的k  
for i=2^(kmax-1)+1:h-2^(kmax-1)  
    for j=2^(kmax-1)+1:w-2^(kmax-1)  
        [maxEh,p]=max(Eh(i,j,:));  
        [maxEv,q]=max(Ev(i,j,:));  
        if maxEh>maxEv  
            maxkk=p;  
        else  
            maxkk=q;  
        end  
        Sbest(i,j)=2^maxkk; %每个像素点的最优窗口大小为2^maxkk  
    end  
end  
%所有Sbest的均值作为整幅图片的粗糙度  
Fcrs=mean2(Sbest);  
end 


%% 第二个指标 Contrast，对比度
%注意这个函数因为涉及到方差，要求输入类型为double，因此我这里在源代码上做了适当的修改  
function Fcon=contrast(graypic) %graypic为待处理的灰度图片  
graypic=double(graypic);%这一句我自己做了修改，否则原博文中的代码不能直接运行  
x=graypic(:); %二维向量一维化  
M4=mean((x-mean(x)).^4); %四阶矩  
delta2=var(x,1); %方差  
alfa4=M4/(delta2^2); %峰度  
delta=std(x,1); %标准差  
Fcon=delta/(alfa4^(1/4)); %对比度  
end  

%% 第三个指标 Directionality，方向度
%sita为各像素点的角度矩阵，在线性度中会用到，所以这里作为结果返回  
function [Fdir,sita]=directionality(graypic)  
[h w]=size(graypic);  
%两个方向的卷积矩阵  
GradientH=[-1 0 1;-1 0 1;-1 0 1];  
GradientV=[ 1 1 1;0 0 0;-1 -1 -1];  
%卷积，取有效结果矩阵  
MHconv=conv2(graypic,GradientH);  
MH=MHconv(3:h,3:w);  
MVconv=conv2(graypic,GradientV);  
MV=MVconv(3:h,3:w);  
%向量模  
MG=(abs(MH)+abs(MV))./2;  
%有效矩阵大小  
validH=h-2;  
validW=w-2;  
%各像素点的方向  
for i=1:validH  
    for j=1:validW  
        sita(i,j)=atan(MV(i,j)/MH(i,j))+(pi/2);  
    end  
end  
n=16;  
t=12;  
Nsita=zeros(1,n);  
%构造方向的统计直方图  
for i=1:validH  
    for j=1:validW  
        for k=1:n  
            if sita(i,j)>=(2*(k-1)*pi/2/n) && sita(i,j)<((2*(k-1)+1)*pi/2/n) && MG(i,j)>=t  
                Nsita(k)=Nsita(k)+1;  
            end  
        end  
    end  
end  
for k=1:n  
    HD(k)=Nsita(k)/sum(Nsita(:));  
end  
%假设每幅图片只有一个方向峰值，为计算方便简化了原著  
[maxvalue,FIp]=max(HD);  
Fdir=0;  
for k=1:n  
    Fdir=Fdir+(k-FIp)^2*HD(k);%公式与原著有改动  
end  
end  

%% 第四个指标 Linelikeness，线性度
%image=rgb2gray(imread('example.jpg'));  
%Flin=linelikeness(image,sita,4) %sita为directionality.m返回的结果  
function Flin=linelikeness(graypic,sita,d) %d为共生矩阵计算时的像素间隔距离  
n=16;  
[h,w]=size(graypic);  
%构造方向共生矩阵  
PDd1=zeros(n,n);  
PDd2=zeros(n,n);  
PDd3=zeros(n,n);  
PDd4=zeros(n,n);  
PDd5=zeros(n,n);  
PDd6=zeros(n,n);  
PDd7=zeros(n,n);  
PDd8=zeros(n,n);  
for i=d+1:h-d-2  
    for j=d+1:w-d-2  
        for m1=1:n  
            for m2=1:n  
                %下方向   
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i+d,j)>=(2*(m2-1)*pi/2/n) && sita(i+d,j)<((2*(m2-1)+1)*pi/2/n))  
                    PDd1(m1,m2)=PDd1(m1,m2)+1;  
                end  
                %上方向  
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i-d,j)>=(2*(m2-1)*pi/2/n) && sita(i-d,j)<((2*(m2-1)+1)*pi/2/n))  
                    PDd2(m1,m2)=PDd2(m1,m2)+1;  
                end  
                %右方向  
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i,j+d)>=(2*(m2-1)*pi/2/n) && sita(i,j+d)<((2*(m2-1)+1)*pi/2/n))  
                    PDd3(m1,m2)=PDd3(m1,m2)+1;  
                end  
                %左方向  
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i,j-d)>=(2*(m2-1)*pi/2/n) && sita(i,j-d)<((2*(m2-1)+1)*pi/2/n))  
                    PDd4(m1,m2)=PDd4(m1,m2)+1;  
                end  
                %右下方向  
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i+d,j+d)>=(2*(m2-1)*pi/2/n) && sita(i+d,j+d)<((2*(m2-1)+1)*pi/2/n))  
                    PDd5(m1,m2)=PDd5(m1,m2)+1;  
                end  
                %右上方向  
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i-d,j+d)>=(2*(m2-1)*pi/2/n) && sita(i-d,j+d)<((2*(m2-1)+1)*pi/2/n))  
                    PDd6(m1,m2)=PDd6(m1,m2)+1;  
                end  
                %左下方向  
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i+d,j-d)>=(2*(m2-1)*pi/2/n) && sita(i+d,j-d)<((2*(m2-1)+1)*pi/2/n))  
                    PDd7(m1,m2)=PDd7(m1,m2)+1;  
                end  
                %左上方向  
                if (sita(i,j)>=(2*(m1-1)*pi/2/n) && sita(i,j)<((2*(m1-1)+1)*pi/2/n)) && (sita(i-d,j-d)>=(2*(m2-1)*pi/2/n) && sita(i-d,j-d)<((2*(m2-1)+1)*pi/2/n))  
                    PDd8(m1,m2)=PDd8(m1,m2)+1;  
                end  
            end  
        end  
    end  
end  
f=zeros(1,8);  
g=zeros(1,8);  
for i=1:n  
    for j=1:n  
        f(1)=f(1)+PDd1(i,j)*cos((i-j)*2*pi/n);  
        g(1)=g(1)+PDd1(i,j);  
        f(2)=f(2)+PDd2(i,j)*cos((i-j)*2*pi/n);  
        g(2)=g(2)+PDd2(i,j);  
        f(3)=f(3)+PDd3(i,j)*cos((i-j)*2*pi/n);  
        g(3)=g(3)+PDd3(i,j);  
        f(4)=f(4)+PDd4(i,j)*cos((i-j)*2*pi/n);  
        g(4)=g(4)+PDd4(i,j);  
        f(5)=f(5)+PDd5(i,j)*cos((i-j)*2*pi/n);  
        g(5)=g(5)+PDd5(i,j);  
        f(6)=f(6)+PDd6(i,j)*cos((i-j)*2*pi/n);  
        g(6)=g(6)+PDd6(i,j);  
        f(7)=f(7)+PDd7(i,j)*cos((i-j)*2*pi/n);  
        g(7)=g(7)+PDd7(i,j);  
        f(8)=f(8)+PDd8(i,j)*cos((i-j)*2*pi/n);  
        g(8)=g(4)+PDd8(i,j);  
    end  
end  
tempM=f./g;  
Flin=max(tempM);%取8个方向的线性度最大值作为图片的线性度  
end  


%% 第五个指标 Regularity，规则度  
%image=rgb2gray(imread('example.jpg'));  
%Freg=regularity(image,64)  
function Freg=regularity(graypic,windowsize) %windowsize为计算规则度的子窗口大小  
[h,w]=size(graypic);  
k=0;  
for i=1:windowsize:h-windowsize  
    for j=1:windowsize:w-windowsize  
        k=k+1;  
        crs(k)=coarseness(graypic(i:i+windowsize-1,j:j+windowsize-1),4); %粗糙度  
        con(k)=contrast(graypic(i:i+windowsize-1,j:j+windowsize-1)); %对比度  
        [dire(k),sita]=directionality(graypic(i:i+windowsize-1,j:j+windowsize-1));%方向度  
        lin=linelikeness(graypic(i:i+windowsize-1,j:j+windowsize-1),sita,4)*10; %线性度，*10与crs、con、dire同量级化  
    end  
end 
Dcrs=std(crs,1);  
Dcon=std(con,1);  
Ddir=std(dire,1);  
Dlin=std(lin,1);
Freg=1-(Dcrs+Dcon+Ddir+Dlin)/4/100;
end