getimagecodewords.m

特征向量 主要在MATLAB下实现.用于图像的PCA之前的准备工作

function [CodeWordsNorm,CodeWordsOrg,Blocks,MM,NN,Indx,DD,PCA,BlockCont]=GetImageCodeWords(F,K,SVth,Num,CIEflag)
%---------------------------------------------
%   [CodeWordsNorm,CodeWordsOrg,Blocks,MM,NN,Indx,DD,PCA,BlockCont]=GetImageCodeWords(F,K,SVth,Num,CIEflag)  
%       is called by ExamineCodeWords and CBIRbookGenerator. 
%       It returns codewords derived from features of the input image F,  
%       which are the most significant 6 features derived from features 
%       of image blocks. The features include the 1st to the 4th moments,
%       Tamura coarseness, modified Tamura contrast, 4 color features based on
%       the CIE L*a*b* color representation, and sigular values of the image
%       blocks. Blocks without enough information are discarded. This is
%       achieved by thresholding the block's SV1/SV2. When SV1/SV2 is greater
%       than, say 20, the block is usually too flat to be meaningful in CBIR. 
%       The basic 12 features are reduced to 6 using PCA. The reserved PCA are
%       normalized to the standard deviation within each image. The first column
%       corresponding to the most important components in the PCA domain.
%
%   F:      input image
%   K:      size of block (default=16)
%   SVth:   threshold of SV1/SV2 above which the block is not used (deault=20)
%   Num:    No. of principal components reserved in the output (default=6)
%   CIEflag: 1 for CIE L*a*b* (default), 0 for YCbCR 
%   CodeWordsNorm: T1-by-6 matrix containing T1 code words of 6 numbers,
%           containing normalized codewords (with respect to std). 
%           T1 is th number of informative blocks
%   CodeWordsOrg: T1-by-6 matrix containing T1 code words of 6 numbers,
%           containing original codewords before normalization
%   Blocks is a KxKxLLxT array containing T blocks of the image. LL=3 for color, L=1 for B&W
%   MM,NN:  size of the trimmed image (integer multiples of K)
%   Indx:   positions of the useful T1 blocks (in Blocks) whose SV1/SV2 < SVth
%   DD:     eigenvalues of covariance matrix (in PCA)
%   PCA:    PCA of feature vectors (empty element removed)
%   BlockCont: contrast calculated in the block as a whole
%
%   Copyright(c): S. Wang, January 2007
%------------------------------------------------
if     nargin==1, CIEflag=1;Num=6;SVth=20;K=16;
elseif nargin==2, CIEflag=1;Num=6;SVth=20; 
elseif nargin==3, CIEflag=1;Num=6;
elseif nargin==4, CIEflag=1; 
end

[M,N,LL]=size(F);
G=double(F);

if LL==3, F1=rgb2ycbcr(F); else F1=F; end
B=double(F1(:,:,1));        % if color, take the luminance
Lm=fix(M/K);
Ln=fix(N/K);
Total=Lm*Ln;                % total number of blocks
MM=Lm*K;NN=Ln*K;

F=F(1:MM,1:NN,:);           % image trimmed to integer multiples of blocks
G=G(1:MM,1:NN,:);           % double version of F

CoarsenessScale=32;
[Coarse,Contrast]=TamuraImageCoarseness(F,CoarsenessScale);
Win=9;
msk=ones(Win,Win)/(Win*Win);
Contrast=conv2(Contrast,msk,'s');
TH=20;Q=10;[HD,GradAng]=TamuraDirectionality(F,TH,Q,0);   % Global directionality
if CIEflag,
    CIE=rgb2lab(F);QQ=10;
else
    if LL==3
        CIE=double(rgb2ycbcr(F));
    else
        CIE=zeros(MM,NN,3);
    end    
    QQ=1;
end
CIEa=round(CIE(:,:,2)/QQ);
CIEb=round(CIE(:,:,3)/QQ);

Features=zeros(Total,12);       % Features can hold 12 features for each image block
BlockCont=zeros(Total,2);

Blocks=zeros(K,K,LL,Total);     % hold image Lm*Ln blocks
b=zeros(K,K,Total);             % hold Lm*Ln luminance blocks
c=zeros(K,K,Total);             % hold Lm*Ln coarseness blocks
d=zeros(K,K,Total);             % hold Lm*Ln directionality blocks
cont=zeros(K,K,Total);          % hold Lm*Ln contrast blocks
ciea=zeros(K,K,Total);          % hold Lm*Ln quantized CIE a blocks
cieb=zeros(K,K,Total);          % hold Lm*Ln quantized CIE b blocks

for m=1:Lm
    for n=1:Ln
        t=n+(m-1)*Ln;
        row=(1:K)+(m-1)*K;
        col=(1:K)+(n-1)*K;
        Blocks(:,:,:,t)=G(row,col,:);       % Take a block of image
        b(:,:,t)=B(row,col);                % Take a luminance block
        c(:,:,t)=Coarse(row,col);           % Take a Coarseness block
        d(:,:,t)=round(GradAng(row,col));   % Take average angle block
        cont(:,:,t)=Contrast(row,col);      % Take a contrast block
        ciea(:,:,t)=CIEa(row,col);          % Take a quantized CIE a* block
        cieb(:,:,t)=CIEb(row,col);          % Take a quantized CIE b* block
    end
end

for t=1:Total
    % Features(:,1) - ratio between the first 2 SVs: SV(1)/SV(2)
    [dummy,S,dummy]=svd(b(:,:,t));Features(t,1)=min(round(S(1,1)/(S(2,2)+eps)),SVth);
    
    if Features(t,1)<SVth        
        % Features(:,2) - Average brightness
        bb=b(:,:,t);Features(t,2)=round(mean(bb(:)));
        
        % Features(:,3) - Standard deviation of brightness
        bb=b(:,:,t);Features(t,3)=round(std(bb(:)));
        
        % Features(:,4) - skewness times 10
        bb=b(:,:,t);
        Features(t,4)=round(skewness(bb(:))*10);
        
        % Features(:,5) - kurtosis
        Features(t,5)=round(kurtosis(bb(:)));

        % Features(:,6) - mean Tamura coarseness
        cc=c(:,:,t);
        Features(t,6)=round(mean(cc(:)));
        
        % Features(:,7) - Tamura contrast
        contr=cont(:,:,t);
        BlockCont(t,1)=round(mean(contr(:)));           % contrast derived from coarseness
        BlockCont(t,2)=TamuraBlockContrast(b(:,:,t));   % originally proposed Tamura contrast
        Features(t,7)=BlockCont(t,2);
        
        % Features(:,8) - dominant direction in degrees
        [dummy,dummy,dd]=find(d(:,:,t)); % remove zeros
        [HH,XX]=hist(dd,180); % HH is the histogram. XX gives bin position
        [HHdominant,AngDominant]=dominant(HH,XX);
        Features(t,8)=round(AngDominant);
        
        % Features(:,9) - average of CIE a components
        AA=ciea(:,:,t);
        Features(t,9)=round(mean(AA(:)));
        
        % Features(:,10) - variance of CIE a components
        Features(t,10)=round(var(AA(:)));
        
        % Features(:,11) - average of CIE b components
        BB=cieb(:,:,t);
        Features(t,11)=round(mean(BB(:)));
        
        % Features(:,12) - variance of CIE b components
        Features(t,12)=round(var(BB(:)));
    end
end

ft=zeros(Total,12);
ndx=zeros(Total,1);
k=0;
for t=1:Total
    if Features(t,1)<SVth
        k=k+1;
        ndx(k,1)=t;
        ft(k,:)=Features(t,:);
    end
end
Indx=ndx(1:k,1);

feat=ft(1:k,:);                         % Remove "empty blocks"
Sigma=cov(feat);
[A,D]=eig(Sigma);
DD=diag(D);		                        % A is KLT matrix. DD contains eigenvalues in
PCA=feat*A;                             % PCA of feat

CodeWordsOrg=PCA(:,12:-1:12-Num+1);     % Take the Num most important components

st=std(CodeWordsOrg);
ST=zeros(size(CodeWordsOrg));
for t=1:length(CodeWordsOrg)
    ST(t,:)=st;
end
CodeWordsNorm=round(100*CodeWordsOrg./(ST+eps));  % normalize to the std

%=========================================================================
function [D,XX]=dominant(y,x,TH)
%---------------------------------------------------
%   [D,XX]=dominant(y,x,TH) returns the dominant value in the input matrix y
%       corresponding to its abscissa x, if the maximum value in y is greater 
%       than the average by a factor of TH. If there is no dominant value, D=-1;
%       Default of TH is 2;
%       XX is position (abscissa) of the dominant component in y
%----------------------------------------------------
if nargin<3, TH=2; end
[dummy,dummy,yy]=find(y);
av=sum(yy(:)/(length(yy(:))+eps));
[D,J]=max(y(:));
XX=x(J);
if D<av*TH, D=-1; end

%==========================================================================
function CIE=rgb2lab(f)
%----------------------------------------------------
%   CIE=rgb2lab(f) converts RGB image f into the CIE L*a*b* color space
%       representation.
%       Ref: Y. Rubner, et al., Empirical Evaluation of Dissimilarity
%       Measures for Color and Texture, Computer Vision and Image
%       Understanding, 84, 2001:25-43
%----------------------------------------------------
T=[0.412453 0.357580 0.180423;0.212671 0.715160 0.072169;0.019334 0.119193 0.950227];
Xn=0.95045;
Yn=1;
Zn=1.088754;
third=1/3;
c=16/116;

f=double(f);
[M,N,L]=size(f);
if L==1
    disp('Warning: This is a monochrome image.');
    g=f;f=zeros(M,N,3);
    f(:,:,1)=g;f(:,:,2)=g;f(:,:,3)=g;
end
CIE=zeros(size(f));
wb=waitbar(0,'RGB to CIE Lab conversion ...');
for m=1:M
    for n=1:N 
        XYZ=T*[f(m,n,1);f(m,n,2);f(m,n,3)];
        if XYZ(2)>0.008856
            CIE(m,n,1)=116*XYZ(2)^third-16;
        else
            CIE(m,n,1)=903.3*XYZ(2);
        end
        CIE(m,n,2)=500*(func(XYZ(1)/Xn,XYZ(2),c)-func(XYZ(2),XYZ(2),c));
        CIE(m,n,3)=200*(func(XYZ(2)/Yn,XYZ(2),c)-func(XYZ(3)/Zn,XYZ(2),c));
    end
    waitbar(m/M);   
end
close(wb);
if L==1
    CIE(:,:,2)=0;CIE(:,:,3)=0;
end

%======================================
function y=func(x,th,c)
if th>0.008856
    y=x^(1/3);
else
    y=7.787*x+c;
end