Matlab特征提取之灰度游程（行程）矩阵-GLRLM

Matlab特征提取之灰度游程（行程）矩阵-GLRLM

前几天参加数学建模比赛，发现全网图形识别检测领域与灰度游程矩阵特征提取的matlab的参考材料很少，下面给出了求灰度游程矩阵以及部分特征提取的代码，供大家交流。

灰度游程（行程）矩阵

一幅图像的灰度游程矩阵反映了图像灰度关于方向，相邻间隔和变化幅度等
综合信息。灰度游程矩阵是对分析影像的局部模式和其排列规则基础之一。灰
度游程矩阵可以实现对一幅图像中同一方向同一灰度值连续出现个数的统计。在
一幅图像上，在某一方向上连续的像素点具有相同的灰度值，灰度游程矩阵就是
通过对这些像素点的分布进行统计得到纹理特征。本文中，我们使用P(i,j|θ)代表
灰度游程矩阵，(i，j)点代表灰度为i的像素在图像θ方向连续出现j次的计数。此
P(i,j|θ)代表灰度游程矩阵，θ一般有0度、45度、90度和135度。
${\mathrm{N}}_{\mathrm{g}}$表示图像上的灰度级数目
$N_{\mathrm{r}}$表示图像上不同游程的数目
${\mathrm{N}}_{\mathrm{p}}$表示图像上像素点的数目：

举个例子：矩阵A为：

$\left(\begin{array}{cccccccc}0& 0& 1& 1& 2& 2& 3& 3\\ 1& 1& 2& 2& 3& 3& 0& 0\\ 2& 2& 3& 3& 0& 0& 1& 1\\ 3& 3& 0& 0& 1& 1& 2& 2\end{array}\right)$

那么P(i,j|0⁰)（θ=0⁰时），水平统计各像素出现的次数：

P(i,j|0)=$\left(\begin{array}{cccc}0& 4& 0& 0\\ 0& 4& 0& 0\\ 0& 4& 0& 0\\ 0& 4& 0& 0\end{array}\right)$

0出现1次的次数为0,0出现两次的次数为4,0出现3次次数为0，以此类推…

同理θ=45⁰：

P(i,j|45⁰)=$\left(\begin{array}{cccc}4& 2& 0& 0\\ 4& 2& 0& 0\\ 4& 2& 0& 0\\ 2& 3& 0& 0\end{array}\right)$

同理θ=90⁰：

P(i,j|90⁰)=$\left(\begin{array}{cccc}8& 0& 0& 0\\ 8& 0& 0& 0\\ 8& 0& 0& 0\\ 8& 0& 0& 0\end{array}\right)$

代码：

// An highlighted block
%-----------------------------------
% author：Courage
% name：tdcup.model2.灰度游程矩阵
%-----------------------------------
function [glrlms,si]= grayrlmatrix(varargin)

[i, offset, nl, gl] = parseinputs(varargin{:});

if gl(2) == gl(1)
    si = ones(size(i));
else
    slope = (nl-1) / (gl(2) - gl(1));
    intercept = 1 - (slope*(gl(1)));
    si = round(imlincomb(slope,i,intercept,'double'));
end

si(si > nl) = nl;
si(si < 1) = 1;
numoffsets = size(offset,1);

if nl ~= 0
   
    for k = 1 : numoffsets
        glrlms{k} = computeglrlm(si,offset(k),nl);
    end
else
    glrlms = [];
end

% --------------------------------------------------------------------
function oneglrlm = computeglrlm(si,offset,nl)

switch offset
    case 1
        % 0 degree
        oneglrlm = rle_0(si,nl);
    case 2
        % 45 degree
        seq = zigzag(si);
        oneglrlm  = rle_45(seq,nl);
    case 3
        % 90 degree
        oneglrlm = rle_0(si',nl);
    case 4
       
        seq = zigzag(fliplr(si));
        oneglrlm = rle_45(seq,nl);
    otherwise
        error('only 4 directions supported')
end


% --------------------------------------------------------------------
function [i, offset, nl, gl] = parseinputs(varargin)

iptchecknargin(1,7,nargin,mfilename);

i = varargin{1};
iptcheckinput(i,{'logical','numeric'},{'2d','real','nonsparse'}, ...
    mfilename,'i',1);

offset = [1;2;3;4];

if islogical(i)
    nl = 2;
else
    nl = 8;
end
gl = getrangefromclass(i);


if nargin ~= 1

    paramstrings = {'offset','numlevels','graylimits'};

    for k = 2:2:nargin

        param = lower(varargin{k});
        inputstr = iptcheckstrs(param, paramstrings, mfilename, 'param', k);
        idx = k + 1;  %advance index to the value portion of the input.
        if idx > nargin
            eid = sprintf('images:%s:missingparametervalue', mfilename);
            msg = sprintf('parameter ''%s'' must be followed by a value.', inputstr);
            error(eid,'%s', msg);
        end

        switch (inputstr)

            case 'offset'

                offset = varargin{idx};
                iptcheckinput(offset,{'logical','numeric'},...
                    {'d','nonempty','integer','real'},...
                    mfilename, 'offset', idx);
                % must be row vector 
                if size(offset,2) ~= 1
                    eid = sprintf('images:%s:invalidoffsetsize',mfilename);
                    msg = 'offset must be an n x 1 array.';
                    error(eid,'%s',msg);
                end
                offset = double(offset);

            case 'numlevels'

                nl = varargin{idx};
                iptcheckinput(nl,{'logical','numeric'},...
                    {'real','integer','nonnegative','nonempty','nonsparse'},...
                    mfilename, 'nl', idx);
                if numel(nl) > 1
                    eid = sprintf('images:%s:invalidnumlevels',mfilename);
                    msg = 'nl cannot contain more than one element.';
                    error(eid,'%s',msg);
                elseif islogical(i) && nl ~= 2
                    eid = sprintf('images:%s:invalidnumlevelsforbinary',mfilename);
                    msg = 'nl must be two for a binary image.';
                    error(eid,'%s',msg);
                end
                nl = double(nl);

            case 'graylimits'

                gl = varargin{idx};
                iptcheckinput(gl,{'logical','numeric'},{'vector','real'},...
                    mfilename, 'gl', idx);
                if isempty(gl)
                    gl = [min(i(:)) max(i(:))];
                elseif numel(gl) ~= 2
                    eid = sprintf('images:%s:invalidgraylimitssize',mfilename);
                    msg = 'gl must be a two-element vector.';
                    error(eid,'%s',msg);
                end
                gl = double(gl);
        end
    end
end

十类游程矩阵相关的纹理统计特征：

1、短游程优势：Short Run Emphasis(SRE)

$\mathrm{S}\mathrm{R}\mathrm{E}=\frac{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}\left[\frac{\rho (\mathrm{i},\mathrm{i}|\theta )}{{\mathrm{j}}^2}\right]}{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}\rho (\mathrm{i},\mathrm{j}|\theta )}$

2、长游程优势：Long Run Emphasis(LRE)

$\mathrm{L}\mathrm{R}\mathrm{E}=\frac{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}{\mathrm{j}}^2\rho (\mathrm{i},\mathrm{j}|\theta )}{{\sum }_{i=1}^{\mathrm{N}}{\sum }_{j=1}^{{\mathrm{N}}_{\mathrm{r}}}\rho (\mathrm{i},\mathrm{j}|\theta )}$

3、灰度不均匀性：Gray Level Non．Uniformity(GLN)

$\mathrm{GLN}=\frac{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\left[{\sum }_{\mathrm{j}=1}^{N_{\mathrm{r}}}\rho (\mathrm{i},j|\theta )\right]}^2}{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}\rho (\mathrm{i},\mathrm{j}|\theta )}$

4、长游程不均匀性：Run Length Non—Uniformity(RLN)

$\mathrm{R}\mathrm{L}\mathrm{N}=\frac{{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}{\left[{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}\rho (\mathrm{i},\mathrm{j}|\theta )\right]}^2}{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}\rho (\mathrm{i},j|\theta )}$

5、游程百分比：Run Percentage(RP)

$RP={\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_r}\frac{\rho (i,j|\theta )}{N_{\rho }}$

6、低灰度级游程优势：Low Gray Level Run Emphasis(LGLRE)

$\mathrm{L}\mathrm{G}\mathrm{L}\mathrm{R}\mathrm{E}=\frac{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}\left[\frac{\rho (\mathrm{i},j|\theta )}{{\mathrm{i}}^2}\right]}{{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}\rho (\mathrm{i},j|\theta )}$

7、高灰度级游程优势：High Gray Level Run Emphasis(HGLRE)

$\mathrm{H}\mathrm{G}\mathrm{L}\mathrm{R}\mathrm{E}=\frac{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}{\mathrm{i}}^2\rho (\mathrm{i},j|\theta )}{{\sum }_{\mathrm{i}=1}^{{\mathrm{N}}_{\mathrm{g}}}{\sum }_{\mathrm{j}=1}^{{\mathrm{N}}_{\mathrm{r}}}\rho (\mathrm{i},j|\theta )}$

8、短游程低灰度级优势：Short Run Low Gray Level Emphasis(SRLGLE)

$SRLGLE=\frac{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_r}\left[\frac{p(i,j|\theta )}{i^2{j}^2}\right]}{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_r}p(i,j|\theta )}$

9、短游程高灰度级优势：Short Run High Gray Level Emphasis(SRHGLE)

$SRHGLE=\frac{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_r}\left[\frac{p(i,j|\theta )i^2}{j^2}\right]}{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_r}p(i,j|\theta )}$

10、长游程低灰度级优势：Long Run Low Gray Level Emphasis(LRLGLE)

$LRLGLE=\frac{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_r}\left[\frac{p(i,j|\theta )j^2}{i^2}\right]}{{\sum }_{i=1}^{N_g}{\sum }_{j=1}^{N_r}p(i,j|\theta )}$

代码：

// An highlighted block

%--------------------------------------------------------------------------
% this program select a roi, qunatize to lower bit level and computing 
% gray level run length matrix and seven texture parameters viz., 
%    1. short run emphasis (sre) 
%    2. long run emphasis(lre)
%    3. gray level non-uniformity (gln)
%    4. run percentage (rp)
%    5. run length non-uniformity (rln)
%    6. low gray level run emphasis (lgre)
%    7. high gray level run emphasis (hgre)
%--------------------------------------------------------------------------
% author: courage
%--------------------------------------------------------------------------
clc, clear, close all
im=imread('file name with path');
figure
imshow(im)
im1=imcrop(im);
im2=im1(1:128,1:128);
im2=double(im2);
[m,n]=size(im2);
% --------- image quantization to 4 bits (16 gray levels)------------------
imax=max(max(im2));
imin=min(min(im2));
newim=im2-imin;
nmax=max(max(newim));
nmin=min(min(newim));
q=round(nmax/16);
[m,n]=size(newim);
quant=0;
for i=1:m
    for j=1:n
        i = newim(i,j);
        for b = 1:16
            if (i>quant)&(i<=quant+q)
                newim(i,j)=b/16;
                quant=quant+q;
            end            
        end
    end
end
newmax=max(max(newim));
newim1=newim/newmax;
newim2=round(newim1*16)+1;
dir=0; 
dist1=1;
if (dir == 1)
    newim2=newim2';
end
mx = max(max(newim2));
mn = min(min(newim2));
gl = (mx-mn)+1;
[p,q] = size(newim2);
n=p*q;
count=1;
c=1;
col=1;
grl(mx,p)=0;
maxcount(p*q)=0;
mc=0;
%---------------------computing gray level run length matrix---------------
for j=1:p
    for k=1:q-dist1
        mc=mc+1;
        g=newim2(j,k);
        f=newim2(j,k+dist1);
        if (g==f)&(g~=0)
            count=count+1;
            c=count;
            col=count;
            maxcount(mc)=count;
        else grl(g,c)=grl(g,c)+1;col=1;
            count=1;
            c=1;
        end
    end
    grl(f,col)=grl(f,col)+1;
    count=1;
    c=1;
end
i=(mx:mn);
m=grl(mn:mx,:);
m1=m';
maxrun=max(max(maxcount));
s=0;
g(gl)=0;
r(q)=0;
for u=1:gl
    for v=1:q
        g(u)=g(u)+m(u,v);
        s=s+m(u,v);
    end
end
for u1=1:q
    for v1=1:gl
        r(u1)=r(u1)+m1(u1,v1);
    end
end
[dim,dim1]=size(g);
sre=0; lre=0; gln=0; rln=0; rp=0; lgre=0; hgre=0;
for h1= 1:maxrun
    sre=sre+(r(h1)/(h1*h1));
    lre=lre+(r(h1)*(h1*h1));
    rln=rln+(r(h1)*r(h1));
    rp=rp+r(h1);
end
sre1=sre/s;
lre1=lre/s;
rln1=rln/s;
rp1=rp/n;
for h2=1:gl
    gln=(gln+g(h2)^2);
    lgre=lgre+(g(h2)/(h2*h2));
    hgre=hgre+(h2*h2)*g(h2);
end
gln1=gln/s;
lgre1=lgre/s;
hgre1=hgre/s;
clc
% ---------------------------display the parameters------------------------
disp(sprintf('%6.4f',sre1))
disp(sprintf('%6.4f',lre1))
disp(sprintf('%6.4f',gln1))
disp(sprintf('%6.4f',rp1))
disp(sprintf('%6.4f',rln1))
disp(sprintf('%6.4f',lgre1))
disp(sprintf('%6.4f',hgre1))