function MGBIM = compute_GBIM(Y, B, nFrames)
%   No-reference Generalized Block Impairment Metric (GBIM), version 1.0
%   Copyright(c) 2011 Matteo Naccari
%   All Rights Reserved.
%
%   The author is with the Instituto de Telecomunicações,
%   Lisbon, Portugal.
%   email: matteo.naccari@lx.it.pt | matteo.naccari@polimi.it
%
%   Permission to use, copy, or modify this software and its documentation
%   for educational and research purposes only and without fee is hereby
%   granted, provided that this copyright notice and the original authors'
%   names appear on all copies and supporting documentation. This program
%   shall not be used, rewritten, or adapted as the basis of a commercial
%   software or hardware product without first obtaining permission of the
%   authors. The authors make no representations about the suitability of
%   this software for any purpose. It is provided "as is" without express
%   or implied warranty.
%
%   The compute_GBIM software computes the perceptual blockiness metric
%   proposed in H. R. Wu and M. Yuen, "A generalized block-edge impairment
%   metric for video coding", IEEE Signal Process. Lett., vol. 4, no. 11,
%   pp. 317-320, Nov. 1997  
%   The metric works in no-reference modality and it may be parametrized
%   with respect to the square image block dimension.
%   Each image is assumed as a 8-bit gray levels content.
%
%   Input parameters:
%       - I                = Luminance component for the image whose GBIM
%                            is being computed
%       - B                = Block size to be considered in the GBIM
%                            computation
%
%   Output parameters:
%       - MGBIM            = Image level metric value

%   Declarations and assignments
% MGBIM=zeros(length);
% Edited by Bc. Ondrej Zach, dep. od Radioelectronics, FEEC BUT, 2013
[R C]=size(Y(:,:,1));
xi = 81;
lambda = log(1 + sqrt(255 - xi))/log(1 + sqrt(xi));
h = waitbar(0,'Please wait...'); %create progress bar
for k=1:nFrames
    I=single(Y(:,:,k));

%   Dimensions check
if rem(R, B) ~= 0
    error('Number of rows (%d) is not a multiple of block size (%d)', R, B);
end
if rem(C, B) ~= 0
    error('Number of columns (%d) is not a multiple of block size (%d)', R, B);
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Horizontal blockiness
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Compute the Dcf matrix as in equation (1)
right_c = B:B:C-B;
left_c  = B+1:B:C-B+1;

image_col_right = I(:, right_c);
image_col_left  = I(:, left_c);

Dcf = image_col_right - image_col_left;

%   Compute the mu(i,j) and the sigma(i,j) terms
mu    = zeros(R, C/B - 1);
sigma = zeros(R, C/B - 1);
for c = 1:C/B - 1
    left_col  = I(:, (c-1)*B + 1:c*B);
    right_col = I(:, c*B + 1:(c+1)*B);
    mu_left   = mean(left_col, 2);
    mu_right  = mean(right_col, 2);    
    mu(:, c) = (mu_left + mu_right)/2;
    
    sigma_left  = std(left_col, 1, 2);
    sigma_right = std(right_col, 1, 2);
    
    sigma(:, c) = (sigma_left + sigma_right)/2;
    
end

%   Compute the intermediate terms for the weight w
index_low = find(mu <= xi);
index_up = find(mu > xi);

%   arg1 = 1 + sqrt(mu)/(1 + sigma). For mu(i, j) <= xi
arg1 = 1 + sqrt(mu)./(1 + sigma);

%   arg2 = 1 + sqrt(255 - mu)/(1 + sigma). For mu(i, j) > xi
arg2 = 1 + sqrt(255 - mu)./(1 + sigma);

%   low_cond = lambda*log(arg1);
low_cond = lambda * log(arg1);

%   up_cond = log(arg2)
up_cond  = log(arg2);

%   Compute the term w as in equation (3)
w = zeros(R, C/B - 1);
w(index_low) = low_cond(index_low);
w(index_up)  = up_cond(index_up);

%   Compute the Mh matrix as in equation (1)
temp = w.*Dcf;
Mh = sqrt(sum(sum(temp.*temp)));

%   Compute the matrix E as in equation (5)
S = zeros(1, B-1);
for n = 1:B-1
    current = I(:, B + n: B :end);
    next    = I(:, B + n + 1: B :end);
    D       = current - next;
    temp    = w.*D;
    S(n)    = sqrt(sum(sum(temp.*temp)));
end

E = mean(S);

%   Compute MhGBIM
MhGBIM = Mh /E;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Vertical blockiness
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   Compute the vertical Dcf as in equation (1)
top_r = B:B:R-B;
bot_r = B+1:B:R-B+1;

image_row_top = I(top_r, :);
image_row_bot = I(bot_r, :);

Dcf = image_row_top - image_row_bot;

%   Compute the mu(i,j) and the sigma(i,j) terms
mu    = zeros(R/B-1, C);
sigma = zeros(R/B-1, C);
for r = 1:R/B - 1
    top_row = I((r-1)*B + 1:r*B, :);
    bot_row = I(r*B + 1:(r+1)*B, :);
    mu_top  = mean(top_row, 1);
    mu_bot  = mean(bot_row, 1);    
    mu(r, :) = (mu_top + mu_bot)/2;
    
    sigma_top = std(top_row, 1, 1);
    sigma_bot = std(bot_row, 1, 1);
    
    sigma(r, :) = (sigma_top + sigma_bot)/2;
    
end

%   Compute the intermediate terms for the weight w
index_low = find(mu <= xi);
index_up  = find(mu > xi);

%   arg1 = 1 + sqrt(mu)/(1 + sigma). For mu(i, j) <= xi
arg1 = 1 + sqrt(mu)./(1 + sigma);

%   arg2 = 1 + sqrt(255 - mu)/(1 + sigma). for mu(i, j) > xi
arg2 = 1 + sqrt(255 - mu)./(1 + sigma);

%   low_cond = lambda*log(arg1);
low_cond = lambda * log(arg1);

%   up_cond = log(arg2)
up_cond  = log(arg2);

%   Compute the weight w as in equation (3)
w = zeros(R/B - 1, C);
w(index_low) = low_cond(index_low);
w(index_up)  = up_cond(index_up);

%   Compute the Mv matrix as in equation (2)
temp = w.*Dcf;
Mv = sqrt(sum(sum(temp.*temp)));

%   Compute the matrix E as in equation (5)
S = zeros(B-1, 1);
for n = 1:B-1
    current = I(B + n: B :end, :);
    next    = I(B + n + 1: B :end, :);
    D       = current - next;
    temp    = w.*D;
    S(n)    = sqrt(sum(sum(temp.*temp)));
end

E = mean(S);

%   Compute the MvGBIM
MvGBIM = Mv /E;

%   Compute the final MGBIM at the frame level
MGBIM(k) = 0.5*MhGBIM + 0.5*MvGBIM;
 waitbar(k / nFrames, h); %update progress bar
end
close(h) % close progress bar

