function varargout = NivGrisLongPlages(varargin)
% Calcul de la matrice de longueurs de plages d'une image en niveaux de gris
% 
% INPUT:
% - image 1xQ, Q>=8
% - nombre de bits de requantification (1 à 8) pour 2 à 256 niveaux de gris
%
% OUTPUT:
% matrice (1,7) de longueurs de plages
%   + SRE  : SHORT RUN EMPHASIS
%   + LRE  : LONG RUN EMPHASIS
%   + GLN  : GRAY LEVEL NON-UNIFORMITY
%   + RP   : RUN PERCENTAGE
%   + RLN  : RUN LENGTH NON-UNIFORMITY
%   + LGRE : LOW GRAY LEVEL RUN EMPHASI
%   + HGRE : HIGH GRAY LEVEL RUN EMPHASIS
%
% d'après Wout Oude Elferink, 
% 13-5-2015 University Of Twente, The Netherlands
% https://fr.mathworks.com/matlabcentral/fileexchange/
%        52640-gray-level-run-lenght-image-statistics/content//glrlm.m
%
% Plugin for OpenIsaac
% version 1
% © 2017-2018 Arthur Berger / Alain Clément - Université d'Angers

%---------------------------------------------------------------------------------------------------

% Copyright 2007 - 2018 Université d'Angers - Author: Alain Clément
%
% This program is free software; you can redistribute it and/or modify it under the terms of the GNU
% General Public License as published by the Free Software Foundation; either version 3 of the
% License, or any later version. This program is distributed in the hope that it will be useful, but
% WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
% PARTICULAR PURPOSE. See the GNU General Public License for more details. You should have received
% a copy of the GNU General Public License along with this program; if not, write to the Free
% Software Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA.

%---------------------------------------------------------------------------------------------------

                      
% PLUGIN PARAMETERS
if (nargin == 1) && strcmp(varargin{1},'-f')
    % input
	varargout{1} = {'IMG','image', ...
                   'DATA','nombre de bits de quantification'};
    % output
	varargout{2} = {'DATA','matrice longueurs de plages'};
	return
end	

% INPUT
ObjIMG = varargin{1};
nbb = IsaacDATA_get(varargin{2},'Val');

%---------------------------------------------------------------------------------------------------

% test de l'image d'entrée
if (IsaacIMG_get(ObjIMG,'BitClass') < 8) || (IsaacIMG_get(ObjIMG,'DimZ') ~= 1)
	error('image incompatible')
end
% requantification image sur nbb bits (valeurs entre 0 et (2^nbb)-1)
if (nbb ~= IsaacIMG_get(ObjIMG,'NbBits'))
	ObjIMG = IsaacIMG_changebits(ObjIMG,nbb);
end	
img = IsaacIMG_img2mat(ObjIMG); 

% to double
img = double(img)+1; 
% Arthur : +1 car aprés l'algorithme ajoute une bordure de 0 (Padarray). 
% (Les niveaux de gris de l'image deviennent entre 1 et 2^nbb).
% initialize glrlm: p(i,j)
% -  with i the amount of bin values (quantization levels)
quantize = 2^nbb;
% -  with j the maximum run length (because yet unknown, assume maximum length
%    of image)
maximgS = max(size(img));
% -  four different orientations are used (0, 45, 90 and 135 degrees)
p0 = zeros(quantize,maximgS); p45 = p0; p90 = p45; p135 = p90;
% add zeros to the borders
img = padarray(img,[1 1]);
% initialize rotation
img45 = imrotate(img,45);
% find the run length for each quantization level
for i = 1:quantize
    % find the pixels corresponding to the quantization level
    BW = int8(img == i); 
    BWr = int8(img45 == i);    
    % find the start and end points of the run length
    G0e   = (BW(2:end-1,2:end-1)  - BW(2:end-1,3:end)) == 1;
    G0s   = (BW(2:end-1,2:end-1)  - BW(2:end-1,1:end-2)) == 1;
    G45e  = (BWr(2:end-1,2:end-1) - BWr(2:end-1,3:end)) == 1;
    G45s  = (BWr(2:end-1,2:end-1) - BWr(2:end-1,1:end-2)) == 1;
    G90e  = (BW(2:end-1,2:end-1)  - BW(3:end,2:end-1)) == 1;
    G90s  = (BW(2:end-1,2:end-1)  - BW(1:end-2,2:end-1)) == 1;
    G135e = (BWr(2:end-1,2:end-1) - BWr(3:end,2:end-1)) == 1;
    G135s = (BWr(2:end-1,2:end-1) - BWr(1:end-2,2:end-1)) == 1;
    % find the indexes
    G0s = G0s'; G0s = find(G0s(:));
    G0e = G0e'; G0e = find(G0e(:));
    G45s = G45s'; G45s = find(G45s(:));
    G45e = G45e'; G45e = find(G45e(:));
    G90s = find(G90s(:));
    G90e = find(G90e(:));
    G135s = find(G135s(:));
    G135e = find(G135e(:));
    % find the lengths
    lengths0 = G0e - G0s + 1;
    lengths45 = G45e - G45s + 1;
    lengths90 = G90e - G90s + 1;
    lengths135 = G135e - G135s + 1;   
    % fill the matrix
    p0(i,:) = hist(lengths0,1:maximgS);
    p45(i,:) = hist(lengths45,1:maximgS);
    p90(i,:) = hist(lengths90,1:maximgS);
    p135(i,:) = hist(lengths135,1:maximgS);    
end
% add all orientations
p = p0 + p45 + p90 + p135;
% calculate the features
totSum = sum(p(:));
% Arthur: Features independants du + 1 sur le niveau de gris: 
SRE = sum(sum(p,1) ./ ((1:maximgS).^2)) / totSum;
LRE = sum(sum(p,1) .* ((1:maximgS).^2)) / totSum;
RLN = sum(sum(p,1) .^2) / totSum;
RP = totSum / ((size(img,1)-2)*(size(img,2)-2)); % sans la bordure à 0
GLN = sum(sum(p,2) .^2) / totSum;
% Arthur : Features influencés par le + 1 mais sinon division par 0 
% L'influence du + 1 reste minime, sur LGRE et HGRE.
% Erreur trouvée dans l'algorithme de Wout Oude Elferink:
% LGRE = sum(sum(p,2) .* ((1:quantize)'.^2)) / totSum;
% HGRE = sum(sum(p,2) .^2) / totSum;
% Correction de LGRE et HGRE d'aprés publication de A.CHU 1990.
LGRE = sum(sum(p,2) ./((1:quantize)'.^2)) / totSum; 
HGRE = sum(sum(p,2) .*((1:quantize)'.^2)) / totSum;

%---------------------------------------------------------------------------------------------------

% OUTPUT
varargout{1} = IsaacDATA('TABN',  ...
    [SRE,LRE,GLN,RP,RLN,LGRE,HGRE],'');
