clear all;

% Define GLCM-parameters.
w = 41;
G = 8;
dx = 1;
dy = 1;

% Read the image to analysis.
img = imread('zebra_1.tif');

% All zebra images have 2^4 = 16 possible intensities. 'zebra_1.tif' stores
% these using an 8 bits intensity range (the other images uses 4 bits as
% they should do). To convert the read 'zebra_1.tif' to use a 4 bits
% intensity range, use e.g.:
img = uint8(floor(double(img)/16)); % only when using 'zebra_1.tif'

% figure, imshow(img, [0 15]) % notice the black-and-white frame in 'zebra_1.tif'

% Make the histogram (approx.) uniform with G grey levels.
img_std = histeq(img,G);
img_std = round(double(img_std) * (G-1) / double(max(img_std(:))));
% round((0:(G-1))*255/(G-1))
% unique(img_std)'
% figure, imshow(img_std, [0 G-1])

% Index matrices.
i = repmat((0:(G-1))', 1, G);
j = repmat( 0:(G-1)  , G, 1);
% could instead have used MESHGRID:
% [j i] = meshgrid(1:G);
% i = i - 1;
% j = j - 1;

% Compute the GLCM features.
w_half = round((w-1)/2);
glcm_var = zeros(size(img_std));
glcm_ctr = zeros(size(img_std));
glcm_ent = zeros(size(img_std));
for m = 1:size(img,1)
    for n = 1:size(img,2)
        % Compute GLCM.
        m_low  = max(m-w_half, 1);
        n_low  = max(n-w_half, 1);
        m_high = min(m+w_half, size(img,1));
        n_high = min(n+w_half, size(img,2));
        img_std_window = img_std(m_low:m_high, n_low:n_high);
        p = task_glcm(img_std_window, dx, dy, G);

        % Compute GLCM's variance feature.
        mu = mean(img_std_window(:));
        glcm_var(m,n) = sum(sum(p .* ((i-mu).^2)));
        % Computed using for-loops:
        % for i = 0:(G-1)
        %     for j = 0:(G-1)
        %         glcm_var(m,n) = glcm_var(m,n) + p(i+1,j+1) * (i-mu)^2;
        %     end
        % end

        % Compute GLCM's contrast feature.
        glcm_ctr(m,n) = sum(sum(p .* ((i - j).^2)));
        % Computed using for-loops:
        % for i = 0:(G-1)
        %     for j = 0:(G-1)
        %         glcm_ctr(m,n) = glcm_ctr(m,n) + p(i+1,j+1) * (i-j)^2;
        %     end
        % end

        % Compute GLCM's entropy feature (with base 2).
        glcm_ent(m,n) = -sum(p(p>0) .* log2(p(p>0)));
    end
end

% c) Display the feature images and the threshold results.
figure, subplot(1,2,1), imshow(glcm_var, []), title('GLCM Variance'), subplot(1,2,2), imshow(glcm_var > (max(glcm_var(:)) * 0.3))
figure, subplot(1,2,1), imshow(glcm_ctr, []), title('GLCM Contrast'), subplot(1,2,2), imshow(glcm_ctr > (max(glcm_ctr(:)) * 0.12))
figure, subplot(1,2,1), imshow(glcm_ent, []), title('GLCM Entropy' ), subplot(1,2,2), imshow(glcm_ent > (max(glcm_ent(:)) * 0.7))

% d) Compare with first order texture measurements.
std_var = stdfilt(img, ones(w)).^2;
figure, subplot(1,2,1), imshow(std_var, []), title('Variance'), subplot(1,2,2), imshow(std_var > (max(std_var(:)) * 0.25))
std_ent = entropyfilt(img, ones(w));
figure, subplot(1,2,1), imshow(std_ent, []), title('Entropy' ), subplot(1,2,2), imshow(std_ent > (max(std_ent(:)) * 0.75))

std_var_of_img_std = stdfilt(img_std, ones(w)).^2;
figure, subplot(1,2,1), imshow(glcm_var, []), title('GLCM Variance'), subplot(1,2,2), imshow(std_var_of_img_std, []), title('Variance')

% e) Use Law's texture masks.
L3 = [ 1  2  1];
E3 = [-1  0  1];
S3 = [-1  2 -1];

L5 = conv2(L3, L3, 'full');
E5 = conv2(L3, E3, 'full');
S5 = conv2(L3, S3, 'full');
W5 = conv2(E3, S3, 'full'); % not mentioned in the lecture notes
R5 = conv2(S3, S3, 'full');

E5E5 = conv2(E5', E5, 'full');

laws_energy = imfilter(double(img), E5E5, 'symmetric', 'conv');
% figure, imshow(laws_energy, [])
laws_mean_abs_feat = imfilter(abs(laws_energy), ones(35), 'symmetric');
figure, subplot(1,2,1), imshow(laws_mean_abs_feat, []), subplot(1,2,2), imshow(laws_mean_abs_feat > (max(laws_mean_abs_feat(:)) * 0.2))
laws_std_feat  = stdfilt(laws_energy, ones(35));
figure, subplot(1,2,1), imshow(laws_std_feat,      []), subplot(1,2,2), imshow(laws_std_feat      > (max(laws_std_feat(:))      * 0.25))