function [theta_opt, IncBin, IncDots] = Demar_Esti(PosLivCell, BDwidth)

RDdots = PosLivCell;

% BDwidth: Bandwidth of the demarcation line

DL = 8;   % Very long distance


theta_range = 0:0.01:pi;

n_theta = length(theta_range);

IncDots = zeros(1, n_theta);  


for m0 = 1:n_theta
    % hold off

    theta = theta_range(m0);
    theta_norm = theta + pi/2;  % normal angle

    TipA = [DL*cos(theta), DL*sin(theta)];
    TipB = -TipA;

    VtxA1 = TipA + (BDwidth/2)*[cos(theta_norm), sin(theta_norm)];
    VtxA2 = TipA - (BDwidth/2)*[cos(theta_norm), sin(theta_norm)];

    VtxB1 = TipB + (BDwidth/2)*[cos(theta_norm), sin(theta_norm)];
    VtxB2 = TipB - (BDwidth/2)*[cos(theta_norm), sin(theta_norm)];

    Coords = [VtxA1; VtxB1; VtxB2; VtxA2];

    cvx_ord = convhull(Coords);

    CoordsCvx = Coords( cvx_ord(1:4),:);

    % scatter(RDdots(:,1), RDdots(:,2),'r.');

    % hold on

    RecShape = polyshape(CoordsCvx(:,1), CoordsCvx(:,2));
    %     plot(RecShape)
    %     xlim([-8 8])
    %     ylim([-8 8])
    %     axis square
    %     pause

    numInc = sum(isinterior(RecShape, RDdots));
    IncDots(m0) = numInc;

end

[~, minOrd] = min(IncDots);

DemarAng = theta_range(minOrd);

theta_opt = DemarAng;

theta_opt_norm = theta_opt + pi/2;  % normal angle

TipAopt = [DL*cos(theta_opt), DL*sin(theta_opt)];
TipBopt = -TipAopt;

VtxA1opt = TipAopt + DL*[cos(theta_opt_norm), sin(theta_opt_norm)];
VtxB1opt = TipBopt + DL*[cos(theta_opt_norm), sin(theta_opt_norm)];


CoordsOpt = [TipAopt; VtxA1opt; VtxB1opt; TipBopt];

cvx_ord = convhull(CoordsOpt);

CoordsCvxOpt = CoordsOpt( cvx_ord(1:4),:);

EstiArea = polyshape(CoordsCvxOpt(:,1), CoordsCvxOpt(:,2));
% plot(EstiArea)


IncBin = isinterior(EstiArea, RDdots);




