function FeaturePlot(Phi,N,alpha,epsilon,xDim,yDim,x,y,xElCount,yElCount,showFeatureSpines,xy00,featureVars,showFeatureIds)

% 24 colours
% colorlist = {[0.7 0.5 1] [0 0 1] [0.3 0.9 0.1] [1 1 0] [1 0.1 0.2] [1 0.5 0.2] [1 0 0] [0 1 0] [0.6 0.8 1] [0.4 0.1 1] [0.7 0.5 1] [0.8 0.4 0.2] [0.7 0.3 0.7] [0.6 0.4 0.5] [0.9 0.7 0] [0 0 1] [0.3 0.9 0.1] [1 1 0] [1 0.1 0.2] [1 0.5 0.2] [1 0 0] [0 1 0] [0.6 0.8 1] [0.4 0.1 1]};
colorlist = {[0 0 1] [0.3 0.9 0.1] [1 1 0] [1 0.1 0.2] [1 0.5 0.2] [1 0 0] [0 1 0] [0.6 0.8 1] [0.4 0.1 1] [0.7 0.5 1] [0.8 0.4 0.2] [0.7 0.3 0.7] [0.6 0.4 0.5] [0.9 0.7 0] [0 0 1] [0.3 0.9 0.1] [1 1 0] [1 0.1 0.2] [1 0.5 0.2] [1 0 0] [0 1 0] [0.6 0.8 1] [0.4 0.1 1] [0.7 0.5 1]};
% % Colour list for 4 feature example in paper
% colorlist = {[0 0 1] [0 1 0] [1 0 0] [1 1 0] [1 0.5 0.2] [1 0 0] [0 1 0] [0.6 0.8 1] [0.4 0.1 1] [0.7 0.5 1] [0.8 0.4 0.2] [0.7 0.3 0.7] [0.6 0.4 0.5] [0.9 0.7 0] [0 0 1] [0.3 0.9 0.1] [1 1 0] [1 0.1 0.2] [1 0.5 0.2] [1 0 0] [0 1 0] [0.6 0.8 1] [0.4 0.1 1] [0.7 0.5 1]};
% % Colour list for 4 feature example in paper
% colorlist = {[1 0 0] [1 1 0]};

for i=1:N 
    plotdata = Heaviside(Phi{i},alpha,xElCount,yElCount,epsilon);
    Q = logical(round(plotdata*100));
    plotdata(~Q) = NaN;
    plotdata = reshape(plotdata,yElCount+1,xElCount+1);
    surf(x, y,plotdata,'FaceColor',colorlist{mod(i,length(colorlist)+1)},'EdgeColor','none','FaceAlpha',0.4);
    hold on
end
axis equal;axis([0 xDim 0 yDim])
xlabel('x-axis (mm)')
ylabel('y-axis (mm)')
set(gca, 'visible', 'on')
set(gca,'FontSize',18)
box on
grid off
view(0,90)
hold off

%%%% Additional Options %%%%
% grid minor

if showFeatureSpines
    % Linear spines
    hold on
    for i=1:N 
        xy = xy00(featureVars*i-featureVars+1: featureVars*i);
        plot([xy(1), xy(3)],[xy(2), xy(4)],'k--o','LineWidth',1.5,'MarkerFaceColor','k');
    end
    hold off
end

if showFeatureIds
    hold on
    for i=1:N 
        x = (xy00(featureVars*(i-1)+3)-xy00(featureVars*(i-1)+1))/2+xy00(featureVars*(i-1)+1);
        y = (xy00(featureVars*(i-1)+4)-xy00(featureVars*(i-1)+2))/2+xy00(featureVars*(i-1)+2);
        text(x,y,sprintf('%d', i),'FontSize',18);
    end
    hold off
end

end