%Use haralickTextureFeatures from 
% https://www.mathworks.com/matlabcentral/fileexchange/58769-haralicktexturefeatures

fSep=filesep;
p = mfilename("fullpath");%get full path of current script
pString=regexp(p,fSep,'split');
mScript=pString{end};%return the name of the current script

bitLevelsMin=4; %must be selcted > gradMax
bitLevelsMax=8; % bit depth lavels 2^bit = Ng
gradMax=8; % maximum gradient per pixel
dxmax=10;% maximum  distnace between pixels for Harlick measure
nPeriods=1;%how may times is the pattern repeated

%generate image gradients
folder1='figures';
folder2='syntheticImg';
path_figs_out=[ '..' fSep,folder1,fSep];
path_imgs_out=[ '..' fSep,folder1,fSep,folder2,fSep];

featureSelect={'Angular Second Moment (Energy)','Contrast','Correlation','Variance',...
    'Inverse Difference Moment (Homogeneity)','Sum Average','Sum Variance',...
'Sum Entropy','Entropy','Difference Variance','Difference Entropy',...
'Information Measure of Correlation I','Information Measure of Correlation II','Maximal Correlation Coefficient'};

for bitLevels=bitLevelsMin:bitLevelsMax
    for grad=1:gradMax
        for d=1:dxmax
            y=0:grad:2^bitLevels-1;
            y2=repmat(y,[1 10]);
            
            img=repmat(y2,[size(y2,2),nPeriods]);
            img=img';
            image{d,grad,bitLevels}=img;%store each image
            
            path_synthetic_img_out=[path_figs_out, 'syntheticImg',fSep,['bitDepth',num2str(bitLevels)],fSep];
            fNameLocal=[path_synthetic_img_out,mScript,'__img_bitDepth',num2str(bitLevels),...
                '_grad',num2str(grad),'_disp',num2str(d)];
            writematrix(img,[fNameLocal,'.csv']);
            I2 = 255*(img - min(img(:))) ./ (max(img(:)) - min(img(:))); %scale values between 0 and 255    
            I2 = cast(I2,'uint8');
            imwrite(I2,[fNameLocal,'.png'])

             glcm=graycomatrix(img, 'Offset',[d 0],'NumLevels',2.^bitLevels,...
                'GrayLimits',[],'Symmetric',true);
            A(:,d,grad,bitLevels) = haralickTextureFeatures(glcm);
        end
    end 
end 

bitLevels=8;
displ=1:dxmax;
feature=6;
subplot(2,2,1);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,:,grad,bitLevels));
    plot(1:dxmax,xx);
    hold on
    yf6=f6(2.^bitLevels,displ,grad);
    plot(displ,yf6);
end
xlabel('displacement vector');
ylabel(featureSelect{feature});
hold off

% 
feature=7;
subplot(2,2,2);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,:,grad,bitLevels));
    plot(1:dxmax,xx);
    hold on
    yf7=f7(2.^bitLevels,displ,grad);
     plot(displ,yf7);
end
xlabel('displacement vector');
ylabel(featureSelect{feature});
hold off

% 
feature=9;
subplot(2,2,3);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,:,grad,bitLevels));
    plot(1:dxmax,xx);
    hold on
    yf9=f9(2.^bitLevels,displ,grad);
    plot(displ,yf9);
end
xlabel('displacement vector');
ylabel(featureSelect{feature});
hold off

% 
feature=10;
subplot(2,2,4);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,:,grad,bitLevels));
    plot(1:dxmax,xx);
    hold on
    yf10=f10(2.^bitLevels,displ,grad);  
     plot(displ,yf10);
end
xlabel('displacement vector');
ylabel(featureSelect{feature});
hold off

%
fNameLocal=[path_figs_out, mScript,'__Haralick_f6_f7_f9_f10_vs_displ_bidDepth_',num2str(bitLevels)];
exportgraphics(gcf,[fNameLocal,'.eps'],'ContentType','image',...
                'BackgroundColor','none','Resolution',300)

%%%%%%%%%%%%%%%%%%%%%%%%%
%For a fixed distance, plot features f6,f7,f9, and f10 vs number of gray
%levels
bitLevels=bitLevelsMin:bitLevelsMax;
displ=1;

feature=6;
subplot(2,2,1);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,displ,grad,:));
    plot(2.^bitLevels,xx(bitLevels));
    hold on
    yf6=f6(2.^bitLevels,displ,grad);
    plot(2.^bitLevels,yf6);
end
xlabel('bit depth');
ylabel(featureSelect{feature});
hold off

% 
feature=7;
subplot(2,2,2);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,displ,grad,:));
    plot(2.^bitLevels,xx(bitLevels));
    hold on
    yf7=f7(2.^bitLevels,displ,grad);
     plot(2.^bitLevels,yf7);
end
xlabel('bit depth');
ylabel(featureSelect{feature});
hold off

% 
feature=9;
subplot(2,2,3);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,displ,grad,:));
    plot(2.^bitLevels,xx(bitLevels));
    hold on
    yf9=f9(2.^bitLevels,displ,grad);
     plot(2.^bitLevels,yf9);
end
xlabel('bit depth');
ylabel(featureSelect{feature});
hold off

% 
feature=10;
subplot(2,2,4);
for grad=1:round(gradMax/4):gradMax
    xx=squeeze(A(feature,displ,grad,:));
    plot(2.^bitLevels,xx(bitLevels));
    hold on
    yf10=f10(2.^bitLevels,displ,grad);
     plot(2.^bitLevels,yf10);
end
xlabel('bit depth');
ylabel(featureSelect{feature});
hold off

%
fNameLocal=[path_figs_out, mScript,'__Haralick_f6_f7_f9_f10_vs_bidDepth_disp=',num2str(displ)];
exportgraphics(gcf,[fNameLocal,'.eps'],'ContentType','image',...
                'BackgroundColor','none','Resolution',300)