function [Area,CI] = EstimateAreaMC(xpoly,ypoly,maxsize,NbPoints,NbEvaluations,Method,VerboseOutput)

switch Method
    case 'Halton'
        RandomPointsTemp = maxsize .* haltonseq(NbPoints* NbEvaluations,2 );
        k = 0;
        for i = 1 : NbEvaluations
            RandomPoints(1 : NbPoints,2 * i - 1  : 2 * i) = RandomPointsTemp(k * NbPoints + 1 : (k + 1) *NbPoints,1 : 2);
            k = k + 1;
        end
    case 'Sobol'
        RandomPointsTemp = maxsize .* wsobol(NbPoints * NbEvaluations + 64,2 );
        RandomPointsTemp(1:64,:)  =  [];
        k = 0;
        for i = 1 : NbEvaluations
            RandomPoints(1 : NbPoints,2 * i - 1  : 2 * i) = RandomPointsTemp(k * NbPoints + 1 : (k + 1) *NbPoints,1 : 2);
            k = k + 1;
        end

    case 'Standard'
        RandomPoints = maxsize .* rand(NbPoints,2 * NbEvaluations);
    otherwise
        error('Invalid Method');
end

for i = 1: NbEvaluations
    IN (:,i) = inpolygon(RandomPoints(:,2* i -1 ),RandomPoints(:,2* i),xpoly,ypoly);
    [Area,StdFinal,CI] = normfit(maxsize.*maxsize.*IN (:,i));
end;




%IN = inpolygon(RandomPoints(:,1),RandomPoints(:,2),xpoly,ypoly);


hold on;
%[Area,dummy,CI] = normfit ((maxsize^2) .* IN);
if (VerboseOutput)
    h = gcf;
    plot(RandomPoints(IN(:,end),end-1),RandomPoints(IN(:,end),end),'g.','LineWidth',1.5);
    plot(RandomPoints(~IN(:,end),end-1),RandomPoints(~IN(:,end),end),'rx','LineWidth',1.5);
    set(h,'WindowStyle','Docked');
    disp(['Area of the Polygon -> ' num2str(polyarea(xpoly,ypoly))]);
 %   disp(['Estimated Area of the Polygon -> ' num2str( (maxsize^2) * sum(IN)./NbTest)]);
    disp(['Estimated Area of the Polygon  -> ' num2str( Area)]);
end