function [P] = generate_pointcloud(b,varargin)
% GENERATE_POINTCLOUD genera una nuvola di punti p interni al bordo b
% Scritto da Pavan Andrea - 01/03/2022
%       [P] = generate_pointcloud(b,varargin)
% INPUT:
%       b: tensore contenente matrici [x,y,n] con i n bordi
% INPUT OPZIONALI:
%       boolInternal: vettore nx1 true se zona interna, false se esterna
%       nnodes: numero punti per cella [default: 1]
%       lmax: lunghezza massima lato cella [default: 1]
%       nextrasteps: numero passi aggiuntivi decomposizione [default: 1]
% OUTPUT:
%       P: matrice nx2 dei punti
pI = inputParser;
addRequired(pI,'b',@ismatrix);
addParameter(pI,'boolInternal',[],@islogical);
addParameter(pI,'nnodes',1,@isnumeric);
addParameter(pI,'lmax',1,@isnumeric);
addParameter(pI,'nextrasteps',1,@isnumeric);
parse(pI,b,varargin{:});
boolInternal = pI.Results.boolInternal;
nnodes = pI.Results.nnodes;
lmax = pI.Results.lmax;
nextrasteps = pI.Results.nextrasteps;

nb = length(b(1,1,:));
if isempty(boolInternal)
    boolInternal = true(nb,1);
end
nodes = [];
for i=1:nb
    nodes = [nodes; b(:,:,i)];
    if abs(b(1,:,i)-b(end,:,i))<2*eps
        nodes(end,:) = [];
    end
end

%quadtree ricorsivo
p = [0];        %vettore parent: p(i)=j indica che cella i è figlia della cella j
x = [0.5*(min(nodes(:,1)+max(nodes(:,1))))];        %coordinata x centro cella i
y = [0.5*(min(nodes(:,2)+max(nodes(:,2))))];        %coordinata y centro cella i
l = [max(nodes(:,1))-min(nodes(:,1))];       %lunghezza lato cella i
n = [length(nodes(:,1))];       %numero punti dentro cella i
while max(n)>nnodes
    nold = n;
    for j=1:length(nold)
        if nold(j)>nnodes || l(j)>lmax
            %decompongo nodo j
            p = [p; j; j; j; j];
            l = [l; l(j)/2; l(j)/2; l(j)/2; l(j)/2];
            x = [x; x(j)+l(j)/4; x(j)-l(j)/4; x(j)-l(j)/4; x(j)+l(j)/4];
            y = [y; y(j)+l(j)/4; y(j)+l(j)/4; y(j)-l(j)/4; y(j)-l(j)/4];
            n = [n; 
                sum((nodes(:,1)>=x(j)).*(nodes(:,2)>=y(j)).*(nodes(:,1)<=x(j)+l(j)/2).*(nodes(:,2)<=y(j)+l(j)/2));
                sum((nodes(:,1)<=x(j)).*(nodes(:,2)>=y(j)).*(nodes(:,1)>=x(j)-l(j)/2).*(nodes(:,2)<=y(j)+l(j)/2));
                sum((nodes(:,1)<=x(j)).*(nodes(:,2)<=y(j)).*(nodes(:,1)>=x(j)-l(j)/2).*(nodes(:,2)>=y(j)-l(j)/2));
                sum((nodes(:,1)>=x(j)).*(nodes(:,2)<=y(j)).*(nodes(:,1)<=x(j)+l(j)/2).*(nodes(:,2)>=y(j)-l(j)/2))];
            l(j) = -1;
            n(j) = -1;
        end
    end
end


%extra-refining ai bordi
if nextrasteps>0
    pl = 1:length(p);       %vettore indici foglie
    for j=2:length(p)
        %per prima cosa si trovano le foglie
        pl(p(j)) = NaN;
    end
    pl = pl(~isnan(pl))';
    for j=1:length(pl)
        if n(pl(j))==nnodes
            %decocomposizioni ulteriori
            xnew = linspace(x(pl(j))-(l(pl(j))/2)+(l(pl(j))/(2^(nextrasteps+1))), x(pl(j))+(l(pl(j))/2)-(l(pl(j))/(2^(nextrasteps+1))), 2^nextrasteps);
            ynew = linspace(y(pl(j))-(l(pl(j))/2)+(l(pl(j))/(2^(nextrasteps+1))), y(pl(j))+(l(pl(j))/2)-(l(pl(j))/(2^(nextrasteps+1))), 2^nextrasteps);
            for i=1:length(ynew)
                p = [p; pl(j)+0*xnew(:)];
                x = [x; xnew(:)];
                y = [y; ynew(i)+0*xnew(:)];
                l = [l; l(pl(j))/(2^nextrasteps)+0*xnew(:)];
                n = [n; 0*xnew(:)];
            end
            l(pl(j)) = -1;
            n(pl(j)) = -1;
        end
    end
end


%bilanciamento
pl = 1:length(p);       %vettore indici foglie
for j=2:length(p)
    %per prima cosa si trovano le foglie
    pl(p(j)) = NaN;
end
pl = pl(~isnan(pl))';
x = x(pl);
y = y(pl);
l = l(pl);
n = n(pl);
for j=1:length(pl)
    %trovo celle vicine

    %bordo N
    idxf1 = abs(y(j)+l(j)/2 -(y-l/2))<=2*eps;
    idxf2 = x(j)>x-l/2;
    idxf3 = x(j)<x+l/2;
    idxN = idxf1.*idxf2.*idxf3;

    %bordo W
    idxf1 = abs(x(j)-l(j)/2 -(x+l/2))<=2*eps;
    idxf2 = y(j)>y-l/2;
    idxf3 = y(j)<y+l/2;
    idxW = idxf1.*idxf2.*idxf3;

    %bordo S
    idxf1 = abs(y(j)-l(j)/2 -(y+l/2))<=2*eps;
    idxf2 = x(j)>x-l/2;
    idxf3 = x(j)<x+l/2;
    idxS = idxf1.*idxf2.*idxf3;

    %bordo E
    idxf1 = abs(x(j)+l(j)/2 -(x-l/2))<=2*eps;
    idxf2 = y(j)>y-l/2;
    idxf3 = y(j)<y+l/2;
    idxE = idxf1.*idxf2.*idxf3;

    idx = or(or(or(idxN,idxW),idxS),idxE);
    idx = find(idx);

    %suddivido celle vicine troppo grandi
    for i=1:length(idx)
        k = idx(i);
        if l(k)>1.5*l(j)
            l = [l; l(k)/2; l(k)/2; l(k)/2; l(k)/2];
            x = [x; x(k)+l(k)/4; x(k)-l(k)/4; x(k)-l(k)/4; x(k)+l(k)/4];
            y = [y; y(k)+l(k)/4; y(k)+l(k)/4; y(k)-l(k)/4; y(k)-l(k)/4];
            n = [n; 0; 0; 0; 0];

            l(k) = -1;
            n(k) = -1;
        end
    end
end
nold = [];
x = x(n>=0);
y = y(n>=0);
l = l(n>=0);


%generazione pointcloud
P = [];
for j=1:length(l)
    %verifico quali dei 5 punti della cella j sono interni al dominio

    %punto centrale
    Ptrial = [x(j) y(j)];
    isinternal = true;
    for i=1:nb
        if check_internal_points(Ptrial,b(:,:,i))~=boolInternal(i)
            isinternal = false;
        end
    end
    if isinternal
        P = [P; Ptrial];
    end

end


end



%% check_internal_points

function [isInside] = check_internal_points(P,Dom)
% CHECK_INTERNAL_POINTS stabilisce se un punto è all'interno di un poligono
% Scritto da Pavan Andrea - 18/02/2022
%       [isInside] = check_internal_points(P,D)
% INPUT:
%       P: punto [x y]
%       Dom: dominio poligonale chiuso [x(:) y(:)]
% OUTPUT:
%       isInside: booleano true se punto interno, false se esterno

D = Dom;
if norm(D(end,:)-D(1,:))>2*eps
    D = [D; D(1,:)];
end

isInside = false;
count = 0;
for i=1:length(D)-1
    if (D(i,1)-P(1))*(D(i+1,1)-P(1))<=0
        yint = interp1([D(i,1) D(i+1,1)],[D(i,2) D(i+1,2)],P(1));
        if P(2)<=yint
            count = count+1;
        end
    end
end

if mod(count,2)==1
    isInside = true;
end

end
