Tamura纹理特征的matlab实现

本文转自：http://blog.csdn.net/jzwong/article/details/51584535

事实上这篇文章并非原创，代码都是别人写的，可是在我的机子上有些地方不能run，我做了一丁点的修改，所以就把文章设置为原创了。另外，最初的参照博客已经看不到了，我看到的已经是别人转的，所以我有必要贴一下别人的链接以示尊重。

第一个链接： http://blog.sina.com.cn/s/blog_59ead5d90100gx1d.html 该文章已被加密，打不开了幸亏有别人已经转载博主的文章，我的这篇文章就是从第二篇文章那里复制加修改得到的，代码整理放在CSDN博客上看起来更加美观，且亲测可直接运行。

第二个链接：http://blog.sina.com.cn/s/blog_5ae7a1de01012r03.html

对应原文为： Textural Features Corresponding to Visual Perception

下面分模块贴代码，代码是我在第二个链接的基础上稍作修改得到的，亲测可直接运行，没有任何问题。

第一个指标 Coarseness，粗糙度

[plain] view plain copy

 

1. %调用举例：  
2. %image=rgb2gray(imread('example.jpg'));  
3. %f=coarseness(image,5)  
4. function Fcrs = coarseness( graypic,kmax )%graphic为待处理的灰度图像，2^kmax为最大窗口  
5. [h,w]=size(graypic); %获取图片大小  
6. A=zeros(h,w,2^kmax); %平均灰度值矩阵A  
7. %计算有效可计算范围内每个点的2^k邻域内的平均灰度值  
8. for i=2^(kmax-1)+1:h-2^(kmax-1)  
9.     for j=2^(kmax-1)+1:w-2^(kmax-1)  
10.         for k=1:kmax  
11.             A(i,j,k)=mean2(graypic(i-2^(k-1):i+2^(k-1)-1,j-2^(k-1):j+2^(k-1)-1));  
12.         end  
13.     end  
14. end  
15. %对每个像素点计算在水平和垂直方向上不重叠窗口之间的Ak差  
16. for i=1+2^(kmax-1):h-2^(kmax-1)  
17.     for j=1+2^(kmax-1):w-2^(kmax-1)  
18.         for k=1:kmax  
19.             Eh(i,j,k)=abs(A(i+2^(k-1),j,k)-A(i-2^(k-1),j));  
20.             Ev(i,j,k)=abs(A(i,j+2^(k-1),k)-A(i,j-2^(k-1)));  
21.         end  
22.     end  
23. end  
24. %对每个像素点计算使E达到最大值的k  
25. for i=2^(kmax-1)+1:h-2^(kmax-1)  
26.     for j=2^(kmax-1)+1:w-2^(kmax-1)  
27.         [maxEh,p]=max(Eh(i,j,:));  
28.         [maxEv,q]=max(Ev(i,j,:));  
29.         if maxEh>maxEv  
30.             maxkk=p;  
31.         else  
32.             maxkk=q;  
33.         end  
34.         Sbest(i,j)=2^maxkk; %每个像素点的最优窗口大小为2^maxkk  
35.     end  
36. end  
37. %所有Sbest的均值作为整幅图片的粗糙度  
38. Fcrs=mean2(Sbest);  
39. end  

第二个指标 Contrast，对比度

[plain] view plain copy

 

1. %调用举例：  
2. %注意这个函数因为涉及到方差，要求输入类型为double，因此我这里在源代码上做了适当的修改  
3. %image=rgb2gray(imread('example.jpg'));  
4. %f=contrast(image)  
5. function Fcon=contrast(graypic) %graypic为待处理的灰度图片  
6. graypic=double(graypic);%这一句我自己做了修改，否则原博文中的代码不能直接运行  
7. x=graypic(:); %二维向量一维化  
8. M4=mean((x-mean(x)).^4); %四阶矩  
9. delta2=var(x,1); %方差  
10. alfa4=M4/(delta2^2); %峰度  
11. delta=std(x,1); %标准差  
12. Fcon=delta/(alfa4^(1/4)); %对比度  
13. end  

第三个指标 Directionality，方向度

[plain] view plain copy

 

1. %调用举例：  
2. %image=rgb2gray(imread('example.jpg'));  
3. %[Fdir,sita]=directionality(image)  
4.   
5. %sita为各像素点的角度矩阵，在线性度中会用到，所以这里作为结果返回  
6. function [Fdir,sita]=directionality(graypic)  
7. [h w]=size(graypic);  
8. %两个方向的卷积矩阵  
9. GradientH=[-1 0 1;-1 0 1;-1 0 1];  
10. GradientV=[ 1 1 1;0 0 0;-1 -1 -1];  
11. %卷积，取有效结果矩阵  
12. MHconv=conv2(graypic,GradientH);  
13. MH=MHconv(3:h,3:w);  
14. MVconv=conv2(graypic,GradientV);  
15. MV=MVconv(3:h,3:w)  
16. %向量模  
17. MG=(abs(MH)+abs(MV))./2;  
18. %有效矩阵大小  
19. validH=h-2;  
20. validW=w-2  
21. %各像素点的方向  
22. for i=1:validH  
23.     for j=1:validW  
24.         sita(i,j)=atan(MV(i,j)/MH(i,j))+(pi/2);  
25.     end  
26. end  
27. n=16;  
28. t=12;  
29. Nsita=zeros(1,n);  
30. %构造方向的统计直方图  
31. for i=1:validH  
32.     for j=1:validW  
33.         for k=1:n  
34.             if sita(i,j)>=(2*(k-1)*pi/2/n) && sita(i,j)<((2*(k-1)+1)*pi/2/n) && MG(i,j)>=t  
35.                 Nsita(k)=Nsita(k)+1;  
36.             end  
37.         end  
38.     end  
39. end  
40. for k=1:n  
41.     HD(k)=Nsita(k)/sum(Nsita(:));  
42. end  
43. %假设每幅图片只有一个方向峰值，为计算方便简化了原著  
44. [maxvalue,FIp]=max(HD);  
45. Fdir=0;  
46. for k=1:n  
47.     Fdir=Fdir+(k-FIp)^2*HD(k);%公式与原著有改动  
48. end  
49. end  

第四个指标 Linelikeness，线性度

[plain] view plain copy

 

1. %调用举例：  
2. %image=rgb2gray(imread('example.jpg'));  
3. %Flin=linelikeness(image,sita,4) %sita为directionality.m返回的结果  
4. function Flin=linelikeness(graypic,sita,d) %d为共生矩阵计算时的像素间隔距离  
5. n=16;  
6. [h,w]=size(graypic);  
7. %构造方向共生矩阵  
8. PDd1=zeros(n,n);  
9. PDd2=zeros(n,n);  
10. PDd3=zeros(n,n);  
11. PDd4=zeros(n,n);  
12. PDd5=zeros(n,n);  
13. PDd6=zeros(n,n);  
14. PDd7=zeros(n,n);  
15. PDd8=zeros(n,n);  
16. for i=d+1:h-d-2  
17.     for j=d+1:w-d-2  
18.         for m1=1:n  
19.             for m2=1:n  
20.                 %下方向   
21.                 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))  
22.                     PDd1(m1,m2)=PDd1(m1,m2)+1;  
23.                 end  
24.                 %上方向  
25.                 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))  
26.                     PDd2(m1,m2)=PDd2(m1,m2)+1;  
27.                 end  
28.                 %右方向  
29.                 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))  
30.                     PDd3(m1,m2)=PDd3(m1,m2)+1;  
31.                 end  
32.                 %左方向  
33.                 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))  
34.                     PDd4(m1,m2)=PDd4(m1,m2)+1;  
35.                 end  
36.                 %右下方向  
37.                 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))  
38.                     PDd5(m1,m2)=PDd5(m1,m2)+1;  
39.                 end  
40.                 %右上方向  
41.                 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))  
42.                     PDd6(m1,m2)=PDd6(m1,m2)+1;  
43.                 end  
44.                 %左下方向  
45.                 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))  
46.                     PDd7(m1,m2)=PDd7(m1,m2)+1;  
47.                 end  
48.                 %左上方向  
49.                 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))  
50.                     PDd8(m1,m2)=PDd8(m1,m2)+1;  
51.                 end  
52.             end  
53.         end  
54.     end  
55. end  
56. f=zeros(1,8);  
57. g=zeros(1,8);  
58. for i=1:n  
59.     for j=1:n  
60.         f(1)=f(1)+PDd1(i,j)*cos((i-j)*2*pi/n);  
61.         g(1)=g(1)+PDd1(i,j);  
62.         f(2)=f(2)+PDd2(i,j)*cos((i-j)*2*pi/n);  
63.         g(2)=g(2)+PDd2(i,j);  
64.         f(3)=f(3)+PDd3(i,j)*cos((i-j)*2*pi/n);  
65.         g(3)=g(3)+PDd3(i,j);  
66.         f(4)=f(4)+PDd4(i,j)*cos((i-j)*2*pi/n);  
67.         g(4)=g(4)+PDd4(i,j);  
68.         f(5)=f(5)+PDd5(i,j)*cos((i-j)*2*pi/n);  
69.         g(5)=g(5)+PDd5(i,j);  
70.         f(6)=f(6)+PDd6(i,j)*cos((i-j)*2*pi/n);  
71.         g(6)=g(6)+PDd6(i,j);  
72.         f(7)=f(7)+PDd7(i,j)*cos((i-j)*2*pi/n);  
73.         g(7)=g(7)+PDd7(i,j);  
74.         f(8)=f(8)+PDd8(i,j)*cos((i-j)*2*pi/n);  
75.         g(8)=g(4)+PDd8(i,j);  
76.     end  
77. end  
78. tempM=f./g;  
79. Flin=max(tempM);%取8个方向的线性度最大值作为图片的线性度  
80. end  

第五个指标 Regularity，规则度

[plain] view plain copy

 

1. %调用举例：  
2. %image=rgb2gray(imread('example.jpg'));  
3. %Freg=regularity(image,64)  
4. function Freg=regularity(graypic,windowsize) %windowsize为计算规则度的子窗口大小  
5. [h,w]=size(graypic);  
6. k=0;  
7. for i=1:windowsize:h-windowsize  
8.     for j=1:windowsize:w-windowsize  
9.         k=k+1;  
10.         crs(k)=coarseness(graypic(i:i+windowsize-1,j:j+windowsize-1),5); %粗糙度  
11.         con(k)=contrast(graypic(i:i+windowsize-1,j:j+windowsize-1)); %对比度  
12.         [dire(k),sita]=directionality(graypic(i:i+windowsize-1,j:j+windowsize-1));%方向度  
13.         lin=linelikeness(graypic(i:i+windowsize-1,j:j+windowsize-1),sita,4)*10; %线性度，*10与crs、con、dire同量级化  
14.     end  
15. end  
16. %求上述各参数的标准差  
17. Dcrs=std(crs,1);  
18. Dcon=std(con,1);  
19. Ddir=std(dire,1);  
20. Dlin=std(lin,1);<pre name="code" class="plain"><span style="font-family: Arial, Helvetica, sans-serif;">%规则度</span>  

Freg=1-(Dcrs+Dcon+Ddir+Dlin)/4/100;end

第六个指标 Roughness，粗略度

[plain] view plain copy

 

1. %粗略度  
2. Frgh=Fcrs+Fcon;  

这篇文章我读的很费劲，感觉很多地方其实并没有多么明白，只是希望拿这几个指标作为图像纹理“好坏”的衡标准，为调参过程提供依据，仅此而已。

有问题可发邮件至 [e-mail] 讨论交流。注：代码非我所写，我只做丁点修改使之可在我的机器上顺利运行。向原作致谢，如有侵权马上删掉！