function [Phi,PhiMax,H,H_features]=FeatureConstruction(xy00,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount)

nodeCount = length(nodeCoords);

Phi=cell(N,1); 
% H_features = zeros((xElCount+1)*(yElCount+1),1);
H_features = zeros(nodeCount,1);
for i=1:N 
    Phi{i}=tPhi(xy00(featureVars*i-featureVars+1: featureVars*i),nodeCoords(:,1),nodeCoords(:,2),p); 
    H_temp=Heaviside(Phi{i},alpha,xElCount,yElCount,epsilon);
    H_features(:,i) = H_temp(:);
end

tempPhiMax = FeatureAggregation(Phi,N,aggregationSelection);
PhiMax=reshape(tempPhiMax,yElCount+1,xElCount+1);

% Heaviside applied to aggregated features
H=Heaviside(PhiMax,alpha,xElCount,yElCount,epsilon);
H=H(:);

end

%Forming Phi_i for each component 
function [tmpPhi]=tPhi(xy,LSgridx,LSgridy,p)
    at_tht = atan2(xy(4)-xy(2),xy(3)-xy(1));
    L = ((xy(3)-xy(1))^2+(xy(4)-xy(2))^2)^0.5;
    xcrd = LSgridx - xy(1)/2 - xy(3)/2;
    ycrd = LSgridy - xy(2)/2 - xy(4)/2;
    tmpPhi = 1 - (2*L^-1*(xcrd*cos(at_tht) + ycrd*sin(at_tht))).^p - (2*(-xcrd*sin(at_tht) + ycrd*cos(at_tht))/xy(5)).^p;
end