function [score] = FQPath(input_image, kernel_sheet)
%
%Non-reference Focus Quality Assessment (NR-FQA) metric for digital 
%pathology
%
%   [score] = FQPath(input_image, params)
%
%   returns the NR-FQA score of the input digital pathological image.
%   Higher score indicates more blurriness in the input image.
%
%
%   Input(s):
%   'input_image'           the grayscale input image for focus quality 
%                           assessment
%                   
%   'params.kernel_sheets'  the pre-designed human visual system-like 
%                           kernel that extracts the sharpness featue from
%                           the input image                  
%
%   Output(s):
%   'score'                 the NR-FQA score 
%
%
%   Copyright (c) 2018 Mahdi S. Hosseini and Yueyang Zhang
%
%   University of Toronto
%   The Edward S. Rogers Sr. Department of,
%   Electrical and Computer Engineering (ECE)
%   Toronto, ON, M5S3G4, Canada
%   Tel: +1 (416) 978 6845
%   email: mahdi.hosseini@mail.utoronto.ca

l = 16;
amp = 1;
alpha = 2;
beta = 1.5;
h_LPF = generalized_Gaussian_for_fitting([-l:l], amp, alpha, beta);
h_LPF = h_LPF/sum(h_LPF);

% do convolution
i_BP_v = imfilter(input_image, h_LPF(:)', 'symmetric', 'conv');
i_BP_v = imfilter(i_BP_v, kernel_sheet, 'symmetric', 'conv');
i_BP_h = imfilter(input_image, h_LPF(:), 'symmetric', 'conv');
i_BP_h = imfilter(i_BP_h, kernel_sheet', 'symmetric', 'conv');
mask = (i_BP_v>0) & (i_BP_h>0);
%mask = iB;
%
v = [abs(i_BP_v(mask)), abs(i_BP_h(mask))];
[pdf, x] = hist(v(:), 50);
cdf = cumsum(pdf)/sum(pdf);
% find sigma approximate
threshold = .95;
min_val = min(cdf);
max_val = max(cdf);
rng_val = max_val - min_val;
indx = cdf < min_val + threshold*rng_val;
sigma_apprx = x(sum(indx))/max(x);
c = (1-tanh(60*(sigma_apprx-.095)))/4 + 0.09;

p_norm = 1/2;
feature_map = (abs(v(:, 1)).^p_norm + abs(v(:, 2)).^p_norm).^(1/p_norm);

%%
number_of_pixels = round(c*numel(feature_map));
feature_map = sort(feature_map(:), 'descend');
feature_map = feature_map(1: number_of_pixels);

%% iterate moments
val = moment(feature_map, 4);
val = abs(val);
val = -log10(val);
if val == (-inf)
    val = 0;
elseif val == inf
    val = 120;
end
score = val;