function [p,t]=AdvancingFront(p,f,Nf)

% FORMULETTE SGUAPPE PER IL CALCOLO DEL CIRCOCENTRO E CIRCORAGGIO DI UN
% TRIANGONO IN 3D DATI I SUOI 3 VERTICI a, b, c
% 
%     |                                                           |
%     | |c-a|^2 (b-a)x(c-a)]x(b-a) + |b-a|^2 (c-a)x[(b-a)x(c-a)   |
%     |                                                           |
% r = -------------------------------------------------------------,
%                          2 | (b-a)x(c-a) |^2
% 
%         |c-a|^2 (b-a)x(c-a)]x(b-a) + |b-a|^2 (c-a)x[(b-a)x(c-a)
% m = a + ---------------------------------------------------------.
%                            2 | (b-a)x(c-a) |^2


alpha=1.15;

NodiSuperficiali=p;
IndiciNodiSuperficiali=1:size(p,1);
FacceSuperficiali=f;

%---------- INIZIALIZZAZIONE DATI DI SUPERFICIE ------------
Spigoli=[f(:,[1 2]); f(:,[1 3]); f(:,[2 3])];
Spigoli=unique(sort(Spigoli,2),'rows');
NumSpigoli=size(Spigoli,1);
NumFacceSuperficiali=size(FacceSuperficiali,1);
Ls=zeros(NumSpigoli,1);
for i=1:NumSpigoli
    Ls(i)=norm([p(Spigoli(i,1),1)-p(Spigoli(i,2),1) p(Spigoli(i,1),2)-p(Spigoli(i,2),2) p(Spigoli(i,1),3)-p(Spigoli(i,2),3)]);
end
Af=zeros(NumFacceSuperficiali,1);
CircoCentri=zeros(NumFacceSuperficiali,3);
CircoRaggi=zeros(NumFacceSuperficiali,1);
F2S=zeros(NumFacceSuperficiali,3); %Per ogni faccia mi dice quali sono i tre spigoli che la compongono
for i=1:NumFacceSuperficiali
    T1=f(i,1); T2=f(i,2); T3=f(i,3); 
    F2S(i,1)=find(sum(abs(ones(NumSpigoli,1)*sort([T1 T2])-Spigoli),2)==0 );
    F2S(i,2)=find(sum(abs(ones(NumSpigoli,1)*sort([T1 T3])-Spigoli),2)==0 );
    F2S(i,3)=find(sum(abs(ones(NumSpigoli,1)*sort([T2 T3])-Spigoli),2)==0 );
    Af(i)=0.5*sqrt(det([p(T1,2) p(T1,3) 1; p(T2,2) p(T2,3) 1; p(T3,2) p(T3,3) 1])^2+det([p(T1,3) p(T1,1) 1; p(T2,3) p(T2,1) 1; p(T3,3) p(T3,1) 1])^2+det([p(T1,1) p(T1,2) 1; p(T2,1) p(T2,2) 1; p(T3,1) p(T3,2) 1])^2);
    [M I]=max(Ls(F2S(i,:)));
    if (sum(Ls(F2S(i,:)).^2)<2*Ls(F2S(i,I))^2) % se il triangolo è ottuso
        CircoRaggi(i)=Ls(F2S(i,I))/2;
        CircoCentri(i,:)=(p(Spigoli(F2S(i,I),1),:)+p(Spigoli(F2S(i,I),2),:))./2;
    else
        v1=[p(T2,1)-p(T1,1); p(T2,2)-p(T1,2); p(T2,3)-p(T1,3)]; %b-a delle formuline
        v2=[p(T3,1)-p(T1,1); p(T3,2)-p(T1,2); p(T3,3)-p(T1,3)]; %c-a delle formuline
        v_cross=cross(v1,v2);                                    %(b-a)x(c-a) delle formuline
        CC=(norm(v2)^2*cross(v_cross,v1)+norm(v1)^2*cross(v2,v_cross))./(2*norm(v_cross)^2);
        CircoRaggi(i)=norm(CC);
        CircoCentri(i,:)=CC'+[p(T1,1) p(T1,2) p(T1,3)];
    end
end
NumNodiSuperficiali=size(NodiSuperficiali,1);
NormaliNodi=zeros(NumNodiSuperficiali,3);
AreaMediaNodo=zeros(NumNodiSuperficiali,1);
rs=zeros(NumNodiSuperficiali,1);
for i=1:NumNodiSuperficiali
    [IndiciFacceIncidenti App]=find(FacceSuperficiali==i); %trova le facce che incidono su quel nodo
    NormaliNodi(i,:)=sum(Nf(IndiciFacceIncidenti,:),1);
    NormaliNodi(i,:)=NormaliNodi(i,:)./norm(NormaliNodi(i,:));
    AreaMediaNodo(i)=mean(Af(IndiciFacceIncidenti));
    [IndiciSpigoliIncidenti App]=find(Spigoli==i); %trova gli spigoli che incidono su quel nodo
    rs(i)=max(Ls(IndiciSpigoliIncidenti));
end
ho=1.24*alpha*sqrt(AreaMediaNodo);
c=ones(NumNodiSuperficiali,1);
Coppie=nchoosek(IndiciNodiSuperficiali,2);
Coppie=Coppie(find( sqrt(sum((p(Coppie(:,1),:)-p(Coppie(:,2),:)).^2,2))<3*rs(Coppie(:,2)) & rs(Coppie(:,2))<rs(Coppie(:,1)) & acos(dot(NormaliNodi(Coppie(:,1),:),NormaliNodi(Coppie(:,2),:),2))>(pi/3)*ones(size(Coppie,1),1)),:);
for i=1:size(Coppie,1)
    %Devo controllare se il nodo Coppie(i,1) è visibile dal nodo Coppie(i,2)
    
    %Estraggo tutte le facce la cui bounding sphere (leggermente + grossa per stare
    %tranquilli) tocca la RETTA che va da Coppie(i,1) a Coppie(i,2)
    %The distance from P2 to the line through P1 in the direction (a,b,c) is
    %sqrt([c(y2-y1)-b(z2-z1)]^2 + [a(z2-z1)-c(x2-x1)]^2 + [b(x2-x1)-a(y2-y1)]^2) / sqrt(a^2+b^2+c^2). 
    %Nel nostro caso: a=p(Coppie(i,2),1)-p(Coppie(i,1),1)
    %                 b=p(Coppie(i,2),2)-p(Coppie(i,1),2)
    %                 c=p(Coppie(i,2),3)-p(Coppie(i,1),3)
    FacceDaControllare=find(((p(Coppie(i,2),2)-p(Coppie(i,1),2))*(CircoCentri(:,1)-p(Coppie(i,1),1)*ones(NumFacceSuperficiali,1))-(p(Coppie(i,2),1)-p(Coppie(i,1),1))*(CircoCentri(:,2)-p(Coppie(i,1),2)*ones(NumFacceSuperficiali,1))).^2 +((p(Coppie(i,2),1)-p(Coppie(i,1),1))*(CircoCentri(:,3)-p(Coppie(i,1),3)*ones(NumFacceSuperficiali,1))-(p(Coppie(i,2),3)-p(Coppie(i,1),3))*(CircoCentri(:,1)-p(Coppie(i,1),1)*ones(NumFacceSuperficiali,1))).^2 + ((p(Coppie(i,2),3)-p(Coppie(i,1),3))*(CircoCentri(:,2)-p(Coppie(i,1),2)*ones(NumFacceSuperficiali,1))-(p(Coppie(i,2),2)-p(Coppie(i,1),2))*(CircoCentri(:,3)-p(Coppie(i,1),3)*ones(NumFacceSuperficiali,1))).^2 <((p(Coppie(i,2),1)-p(Coppie(i,1),1))^2+(p(Coppie(i,2),2)-p(Coppie(i,1),2))^2+(p(Coppie(i,2),3)-p(Coppie(i,1),3))^2)*(1.01)^2*CircoRaggi.^2);
    [FacceIncidenti1 App]=find(f(FacceDaControllare,:)==Coppie(i,1));
    [FacceIncidenti2 App]=find(f(FacceDaControllare,:)==Coppie(i,2));
    FacceResidue=settdiff(setdiff(FacceDaControllare,FacceIncidenti1),FacceIncidenti2);
    
    
    %Le facce intorno al nodo Coppie(i,1) devono avere il nodo Coppie(i,2)
    %che sta dalla parte della normale
    
    
    %In caso affermativo aggiorno c(Coppie(i,1)) e c(Coppie(i,2))
    
end

t=0;