function [originalAreaSum,deltaAreaSum] = AreaDeviationFromInput(H_sub,xy00,featureVars,N,nodeCoords,p,xElCount,yElCount,alpha,epsilon,aggregationSelection,elementNodes,elSize)
% Check shape matching success
    [~,~,H_new,~]=FeatureConstruction(xy00,featureVars,N,nodeCoords,p,alpha,epsilon,aggregationSelection,xElCount,yElCount);

    subdomainArea(:,1:4) = H_sub(elementNodes(:,2:5,1));
    subdomainArea(:,5) = sum(subdomainArea(:,1:4),2)/4*((elSize)^2);
    originalAreaSum = sum(subdomainArea(:,5));

    deltaH = abs(H_new - H_sub);
    deltaArea(:,1:4) = deltaH(elementNodes(:,2:5,1));
    deltaArea(:,5) = sum(deltaArea(:,1:4),2)/4*((elSize)^2);
    deltaAreaSum = sum(deltaArea(:,5)); 

end