% Inputs

% H_features is the implicit representation of features, an array of size(elements x N)
% alpha is the lower limit used for the feature construction (this could be found using min(H_features))
% N is the number of features
% elementNodes is the nodes associated with each element
% elSize is the side length of the element (assumes square elements)

function [featureSize,uniqueSize,intersectingFeatures,intersectingFeatureSize] = FeatureComparison(H_features,alpha,N,elementNodes,elSize)

%Calculating feature element areas
H_area = zeros(length(elementNodes),N);
for i = 1:N
    feature_I = H_features(:,i);
    Area(:,1:4) = feature_I(elementNodes(:,2:5));
    H_area(:,i) = sum(Area(:,1:4),2)/4*((elSize)^2);
end

intersectingFeatures = cell(N,1);
intersectingFeatureSize = cell(N,1);
uniqueSize = zeros(N,1);

loc = logical(H_area>alpha); %Key for feature elements that are not void (size = elements x N)
featureSize = sum(loc,1)';
featureIdList = 1:N;

for i=1:N
    findIntersecting = loc(:,i); %selecting keys for feature i
    %Unique element summation
    uniqueSize(i) = sum(logical(sum(loc(findIntersecting,:),2)==1)); % select rows in loc where findintersecting is active, if the sum of the row is 1 then it is a unqiue element (all non-void elements represented as 1 in loc) 
    %Size of feature overlaps/intersections
	columnFeatureSum = sum(loc(findIntersecting,:),1); % summing the intersecting elements from other features, only where the feature of interest is non-void
    columnFeatureSum(i) = 0; %set the sum of the feature being tested to zero (the previous line calculates the self intersection)
    %Intersecting feature IDs
    featureSelect = logical(columnFeatureSum); % finding the feature ID of features with a non-zero intersection
    [intersectingFeatureSize{i},I] = sort(columnFeatureSum(featureSelect),'descend'); % ranking the intersecting components, most intersecting to least. 
    identifiedIntersectingFeatures = featureIdList(featureSelect); %The feature list is ordered from 1-N so this is unnecessary, included for potential code updates where the input feature list is not sequential from 1-N
    intersectingFeatures{i} = identifiedIntersectingFeatures(I);
end

end