function [ poincare ] = poincare_features( RR )
% 
% Copyright (C) 2017 
% Shreyasi Datta
% Chetanya Puri
% Ayan Mukherjee
% Rohan Banerjee
% Anirban Dutta Choudhury
% Arijit Ukil
% Soma Bandyopadhyay
% Rituraj Singh
% Arpan Pal
% Sundeep Khandelwal
% 
% 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 2
% of the License, or (at your option) 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., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

x = RR(1:end-1);
y = RR(2:end);

% Mean Stepping Increment of Inter Beat Intervals
L = length(RR);
sumval = 0;
for i = 1 : L-2
    sumval = sumval + sqrt((x(i)-y(i))^2 + (x(i+1)-y(i+1))^2);
end
sumval = sumval / (L-2);
stepping = sumval / (mean(RR));


% Dispersion of points around diagonal line in Poincare Plots
cp = (-RR(i)-RR(L)+2*sum(RR(1:L-1)))/(2*(L-1));
term1 = 0;
for i = 1 :  L - 1
    term1 = term1 + (RR(i)-RR(i+1))^2;
end
term1 = term1/(2*(L-1));

term2 = 0;
for i = 1 : L - 1
   term2 = term2 + abs(RR(i) - RR(i+1)); 
end
term2 = term2/(sqrt(2)*(L-1));
term2 = term2^2;

dispersion = sqrt(term1 - term2)/cp;

% Number of clusters in Poincare Plots - try the mclust based clustering
data = [x y];

nobs = size(data,1);
kdist = zeros(nobs, 1);
D = squareform(pdist(data, 'euclidean'));
MinPts = 1;
for i = 1:nobs
    dtmp     = D(i, :);
    dtmp(i)  = [];
    dtmp     = sort(dtmp, 'ascend');
    kdist(i) = dtmp(MinPts);
end
kdist = sort(kdist, 'descend');
[val, loc] = max(kdist);


[clustLabel, varType] = dbscan(data, 1, kdist(loc+1));
clusters = unique(clustLabel);
if ismember(clusters,0)
    numclust = length(unique(clustLabel)) - 1;
else
    numclust = length(unique(clustLabel));
end
%close 
% 
% for i = 1 : length(clustLabel)
%     for j = 1 : length(clusters)
%         if clustLabel(i) == clusters(j) && clustLabel(i) == 0
%             scatter(x(i),y(i),'k');hold on;
%         elseif clustLabel(i) == clusters(j) && clustLabel(i) == 1
%             scatter(x(i),y(i),'b');hold on;
%         elseif clustLabel(i) == clusters(j) && clustLabel(i) == 2
%             scatter(x(i),y(i),'g');hold on;
%         end
%     end
% end
% hold off

% Reducing mistake of missing a QRS Complex (page 7-8)




% Feature List
poincare = [stepping dispersion numclust];

end

