For preparing GLCM feature dataset over several windows and d :
function GLCM(I,q)
% GLCM code
%odd window size is considered and center pixel is considered
%I is the input image
%q is the desired level of quantization
M=size(I,1);%image size MXN(dim=1)
N=size(I,2);%image size MXN(dim=2)
I=I./(256/q);%grey level quantization
winCounter=0;
winCount=0;
winSize=3;
%while winSize^2 <= M*N
% winCount=winCount+1;
%winSize=winSize+2;
%end
Thetan=4;
while winSize^2<=M*N/4
winCounter=winCounter+1;
winSizeTracker(winCounter)=2*winCounter+1;
wc= ceil(winSize/2);%window center
winTracker=0;
winSize,nWin=(M-2*(wc-1))*(N-2*(wc-1))%total no of windows
%featurev=zeros(3,nWin,8);% 3 for 3 channels, nWin no of windows, 8
no of features
%for each window wk and d glcm0==>0 deg glcm45==> 45 deg
%glcm90==> 90 deg glcm135==> 135 deg
%for d=1
%for n number of windows
%for channel red
for di=1:4% max d = as chosen from autocorrelation function plot
if(di<=floor((2*winCounter+1)/2))
for wini=1:M-2*(wc-1)
for winj=1:N-2*(wc-1)
glcm0=zeros(q+1,q+1);
glcm45=zeros(q+1,q+1);
glcm90=zeros(q+1,q+1);
glcm135=zeros(q+1,q+1);
%for each windows
winTracker=(N-2*(wc-1))*(wini-1)+winj;
%------------------------------------------for thetai=1:4
% theta=0 deg
%glcm
for i=wini:wini+winSize-1
for j=winj:winj+winSize-di-1
%dlmwrite('testmat.dat',' -append') ;
glcm0(I(i,j)+1,I(i,j+di)+1)=glcm0(I(i,j)+1,I(i,j+di)+1)+1;
end
end
% symmetric glcm
glcm0=glcm0+(glcm0)';
glcm0=double(glcm0./sum(sum(glcm0)));
%------------------------------------------
[contrast,homog,dissimilarity,entropy,energy,correlation,variance,std_deviation,max_prob,asm]=FeatureCalculator(glcm0,1);
%code khat bala be nazaram javab nemide va be jash az code ir 
%estefade konid
 GLCM = graycomatrix(I,'offset', [0 1;-1 0; -1 1; -1 -1], 'Symmetric', true,'NumLevels',256);
             GLCMSTS= GLCM_Features4(GLCM,0);
             %%%%%%%%%%%%%%%
             if isnan(contrast)
featurev0.contrast=0;
else
featurev0.contrast=contrast;
end
if isnan(homog)
featurev0.homog=0;
else
featurev0.homog=homog;
end
if isnan(dissimilarity)
featurev0.dissimilarity=0;
else
featurev0.dissimilarity=dissimilarity;
end
if isnan(entropy)
featurev0.entropy=0;
else
featurev0.entropy=entropy;
end
if isnan(energy)
featurev0.energy=0;
else
featurev0.energy=energy;
end
if isnan(correlation)
featurev0.correlation=0;
else
featurev0.correlation=correlation;
end
if isnan(variance)
featurev0.variance=0;
else
featurev0.variance=variance;
end
if isnan(std_deviation)
featurev0.std_deviation=0;
else
featurev0.std_deviation=std_deviation;
end
if isnan(max_prob)
featurev0.max_prob=0;
else
featurev0.max_prob=max_prob;
end
if isnan(asm)
    featurev0.asm=0;
else
featurev0.asm=asm;
end
%-------------------------------------------------------------
% theta=45 deg
%glcm
for i=wini+di:wini+winSize-1
for j=winj:winj+winSize-di-1
%dlmwrite('testmat.dat',' -append') ;
glcm45(I(i,j)+1,I(i-di,j+di)+1)=glcm45(I(i,j)+1,I(idi,
j+di)+1)+1;
end
end
% symmetric glcm
glcm45=glcm45+(glcm45)';
glcm45=double(glcm45./sum(sum(glcm45)));
%------------------------------------------
[contrast,homog,dissimilarity,entropy,energy,correlation,variance,std_deviation,max_prob,asm]=FeatureCalculator(glcm45,2);
%_-------------------------------
if isnan(contrast)
featurev45.contrast=0;
else
featurev45.contrast=contrast;
end
if isnan(homog)
featurev45.homog=0;
else
featurev45.homog=homog;
end
if isnan(dissimilarity)
featurev45.dissimilarity=0;
else
featurev45.dissimilarity=dissimilarity;
end
if isnan(entropy)
featurev45.entropy=0;
else
featurev45.entropy=entropy;
end
if isnan(energy)
featurev45.energy=0;
else
featurev45.energy=energy;
end
if isnan(correlation)
featurev45.correlation=0;
else
featurev45.correlation=correlation;
end
if isnan(variance)
featurev45.variance=0;
else
featurev45.variance=variance;
end
if isnan(std_deviation)
featurev45.std_deviation=0;
else
featurev45.std_deviation=std_deviation;
end
if isnan(max_prob)
featurev45.max_prob=0;
else
featurev45.max_prob=max_prob;
end
if isnan(asm)
featurev45.asm=0;
else
featurev45.asm=asm;
end
%------------------------------------------------------------
% theta=90 deg
%glcm
for i=wini+di:wini+winSize-1
for j=winj:winj+winSize-1
%dlmwrite('testmat.dat',' -append') ;
glcm90(I(i,j)+1,I(i-di,j)+1)=glcm90(I(i,j)+1,I(idi,
j)+1)+1;
end
end
% symmetric glcm
glcm90=glcm90+(glcm90)';
glcm90=double(glcm90./sum(sum(glcm90)));
%------------------------------------------
[contrast,homog,dissimilarity,entropy,energy,correlation,variance,std_deviation,max_prob,asm]=FeatureCalculator(glcm90,3);
%_-------------------------------
if isnan(contrast)
featurev90.contrast=0;
else
featurev90.contrast=contrast;
end
if isnan(homog)
featurev90.homog=0;
else
featurev90.homog=homog;
end
if isnan(dissimilarity)
featurev90.dissimilarity=0;
else
featurev90.dissimilarity=dissimilarity;
end
if isnan(entropy)
featurev90.entropy=0;
else
featurev90.entropy=entropy;
end
if isnan(energy)
featurev90.energy=0;
else
featurev90.energy=energy;
end
if isnan(correlation)
featurev90.correlation=0;
else
featurev90.correlation=correlation;
end
if isnan(variance)
featurev90.variance=0;
else
featurev90.variance=variance;
end
if isnan(std_deviation)
featurev90.std_deviation=0;
else
featurev90.std_deviation=std_deviation;
end
if isnan(max_prob)
featurev90.max_prob=0;
else
    featurev90.max_prob=max_prob;
end
if isnan(asm)
featurev90.asm=0;
else
featurev90.asm=asm;
end
%------------------------------------------------------------
% theta=135 deg
%glcm
for i=wini+di:wini+winSize-1
for j=winj+di:winj+winSize-1
%dlmwrite('testmat.dat',' -append') ;
glcm135(I(i,j)+1,I(i-di,jdi)+
1)=glcm135(I(i,j)+1,I(i-di,j-di)+1)+1;
end
end
% symmetric glcm
glcm135=glcm135+(glcm135)';
glcm135=double(glcm135./sum(sum(glcm135)));
%------------------------------------------
[contrast,homog,dissimilarity,entropy,energy,correlation,variance,std_deviation,max_prob,asm]=FeatureCalculator(glcm135,4);
%_-------------------------------
if isnan(contrast)
featurev135.contrast=0;
else
featurev135.contrast=contrast;
end
if isnan(homog)
featurev135.homog=0;
else
featurev135.homog=homog;
end
if isnan(dissimilarity)
featurev135.dissimilarity=0;
else
    featurev135.dissimilarity=dissimilarity;
end
if isnan(entropy)
featurev135.entropy=0;
else
featurev135.entropy=entropy;
end
if isnan(energy)
featurev135.energy=0;
else
featurev135.energy=energy;
end
if isnan(correlation)
featurev135.correlation=0;
else
featurev135.correlation=correlation;
end
if isnan(variance)
featurev135.variance=0;
else
featurev135.variance=variance;
end
if isnan(std_deviation)
featurev135.std_deviation=0;
else
featurev135.std_deviation=std_deviation;
end
if isnan(max_prob)
featurev135.max_prob=0;
else
featurev135.max_prob=max_prob;
end
if isnan(asm)
featurev135.asm=0;
else
featurev135.asm=asm;
end
%calculate average angular feature value
fangularm(winCounter,di,winTracker).contrast=(featurev0.contrast+featurev45.contrast+featurev90.contrast+featurev135.contrast)/4
fangularm(winCounter,di,winTracker).homog=(featurev0.homog+featurev45.homog+featurev90.homog+featurev135.homog)/4
fangularm(winCounter,di,winTracker).dissimilarity=(featurev0.dissimilarity+featurev45.dissimilarity+featurev90.dissimilarity+featurev135.dissimilarity)/4;
fangularm(winCounter,di,winTracker).entropy=(featurev0.entropy+featurev45.entropy+featurev90.entropy+featurev135.entropy)/4;
fangularm(winCounter,di,winTracker).energy=(featurev0.energy+featurev45.energy+featurev90.energy+featurev135.energy)/4;
fangularm(winCounter,di,winTracker).correlation=(featurev0.correlation+featurev45.correlation+featurev90.correlation+featurev135.correlation)/4;
fangularm(winCounter,di,winTracker).variance=(featurev0.variance+featurev45.variance+featurev90.variance+featurev135.variance)/4;
fangularm(winCounter,di,winTracker).std_deviation=(featurev0.std_deviation+featurev45.std_deviation+featurev90.std_deviation+featurev135.std_deviation)/4;
fangularm(winCounter,di,winTracker).max_prob=(featurev0.max_prob+ featurev45.max_prob+ featurev90.max_prob+featurev135.max_prob)/4;
fangularm(winCounter,di,winTracker).asm=(featurev0.asm+featurev45.asm+featurev90.asm+ featurev135.asm)/4;
end
%------------------------------------------------------------
end
end
end
fid=fopen('results101.txt','a+');
fprintf(fid,'For window size %d\n',winSizeTracker(winCounter));
fprintf(fid,'---------------------------------------\n');
for d=1:4
if(d<=floor((2*winCounter+1)/2))
fprintf(fid,'For d= %d\n',d);
fprintf(fid,'Window no|----------Contrast group-------------|--------Orderliness Group--------------------|--DescriptiveStatistics Group--------|\n');
fprintf(fid,' contrast |dissimilarity |homogeniety| asm | energy |max probability| entropy |variance |stddeviation |correlation |\n\n');
for w=1:winTracker
fprintf(fid,'%d\t\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\t%f\n',w,fangularm(winCounter,d,w).contrast,fangularm(winCounter,d,w).dissimilarity,fangularm(winCounter,d,w).homog,fangularm(winCounter,d,w).asm,fangularm(winCounter,d,w).energy,fangularm(winCounter,d,w).max_prob,fangularm(winCounter,d,w).entropy,fangularm(winCounter,d,w).variance,fangularm(winCounter,d,w).std_deviation,fangularm(winCounter,d,w).correlation);
end
end
end
fclose(fid);
winSize=winSize+2;
end